수치해석(4) Roots: Bracketing Methods

 

 

1. How to find Roots

 

(1) Roots

"Roots" 문제는 어떤 함수 f()가 하나 이상의 종속 변수 x의 term으로 쓰여질 때 발생한다. f(x) = 0의 해다. 

 

 

Question : 선형 방정식의 경우, 각 방정식은 Ax = b로 표현 가능한가?

 

(2) Graphical Methods

식 f(x) = 0의 근 추정치를 구하는 간단한 방법은 함수를 그리고, x축과 교차하는 위치를 관찰하는 Graphical Methods다. 함수를 그래프로 그리면 루트가 어디에 있는지, 그리고 루트 찾기 방법 중 일부가 실패한 위치임을 알 수 있다. 일부 예외 상황도 있을 수 있는데, [1] 함수가 x축에 접선일 때 발생하는 다중근, 끝점은 반대 부호이고 짝수 개의 축 절편이 있다. [2] 반대 부호의 끝점이 짝수 개의 근을 괄호로 묶는 불연속 함수가 있다. 

 

 

 

2. Brackeint Methods

 

Bracketing Methods는 루트가 루트의 양쪽에 있는 두 가지의 범위 지정 값을 초기 추측하는 것이다. 즉, 해가 있을 법한 지점의 시작 지점을 x_i로, 끝 지점을 x_u로 두고 이 간격을 점차 좁혀나가면서 해를 찾는다. f(x_i)f(x_u) < 0이다.

또한 incremental search는 균등한 간격으로 함수의 값을 테스트하고, 이웃한 점 사이의 함수 부호 변화를 식별하여 범위가 되어줄 brackets를 찾는다. 증분 길이가 너무 크면, 두 점 내에서 짝수 개의 근을 포착하기 때문에 brackets이 누락될 수 있다. 또한 길이가 너무 작으면, 검색에 시간이 많이 걸릴 수 있다. incremental search에서는 간격에 관계 없이 even-multiplicity roots를 포함하는 brackets를 찾을 수 없다. 

 

 

예제를 풀어보겠다. f(x) = sin(10x) + cos(3x)일 때 incremental search를 수행하자. 

 

#% 분석할 함수, 구간의 최대최소값, 샘플 개수를 입력으로 받는다. 
function xb = incsearch(func, xmin, xmax, ns)
if nargin<3, error('at least 3 arguments required'), end
if nargin<4, ns = 50; end
#% 주어진 구간 [xmin, xmax]에서 'ns'개의 등간격 샘플을 생성한다. 
x=linspace(xmin, xmax, ns);
f=func(x);
#% nb : 브래킷 개수를 나타내는 변수, xb = 브래킷 구간 저장하는 배열. 
nb=0; xb=[]; 
for k=1:length(x)-1
	#% 둘의 부호가 다른 경우, 브래킷이 발견된 것으로 판단. 
    if sign(f(k)) ~= sign(f(k+1))
    	#% nb 값을 증가시키고, 'xb'에 브래킷의 구간을 저장. 
        nb = nb+1;
        xb(nb,1) = x(k);
        xb(nb,2) = x(k+1);
    end
end
#% 브래킷이 없을 경우 메세지를 표시. 
if isempty(xb)
    disp('no brackets found')
    disp('check interval or increase ns')
else
    disp('number of brackets')
    disp(nb)
end

 

 

(1) Bisection Methods

Bisection Methods는 incremental search의 변형 방법으로, 구간을 항상 반으로 나눈다는 특징이 있다. 구간에 걸쳐 부호가 변경되면, 즉 해당 구간 안에 해가 있다는 의미이므로 중간점의 함수값이 평가된다. absolute error는 각 반복에 대해 2배씩 감소한다. 

 

 

Bisection error는 공정 시작 시의 절대 오차와 반복 횟수에만 의존한다. 또한 required number of iterations, 특정 절대 오차를 얻기 위해 필요한 반복 횟수는 초기 추측을 기반으로 계산할 수 있다. 

 

Approximate percent relative error

 

Bisection error : 첫번째는 최대 에러, 두번째는 n번째 iteration 후 최대 에러.

 

Required number of iterations

 

예제를 풀어보겠다. We need to determine the mass at which a bungee jumper’s free-fall velocity exceeds 36 m/s after 4 s of free fall given a drag coefficient of 0.25 kg/m.  

 

#% 이분법을 사용하여 'func'의 근사 해를 찾자. 
#% es : 허용 오차, maxit : 최대 반복 횟수, ea : 추정 상대 오차, iter : 반복 횟수. 
function [root, fx, ea, iter] = bisect(func, xl, xu, es, maxit, varargin)
if nargin<3, error('nonononono !! '), end
#% 초기 구간 끝부분들의 함수 값 부호가 동일하면, 근을 찾을 수 없다. 오류 발생. 
test = func(xl, varargin{:})*func(xu, varargin{:});
if test > 0, error('no sign change'), end
if nargin<4|isempty(es), es=0.001; end
if nargin<5|isempty(maxit), maxit=50; end
iter = 0; xr = xl; ea = 100;
while(1)
    xrold=xr;
    #% xr : 현재 구간의 중점
    xr=(xl+xu)/2;
    iter=iter+1;
    #% 'xr'이 0이 아니라면, 추정 상대 오차 'ea'를 계산. 
    if xr~=0, ea=abs((xr-xrold)/xr)*100;end
    test=func(xl, varargin{:})*func(xr, varargin{:});
    if test<0
        xu=xr;
    elseif test>0
        xl=xr;
    else
        ea=0;
    end
    if ea<=es|iter>=maxit, break, end
end
#% 최종 근사 해 'xr'과 근사 해에서의 함수 값 'fx'를 반환. 
root=xr; fx=func(xr, varargin{:});

 

 

(2) False Position Methods

False Position Methods는 또다른 bracket method로, 브래킷을 반으로 분할하는 것이 아니라 끝점을 직선으로 연결하고 직선(Xr)의 절편 위치를 결정하여 다음 추측을 결정한다. 그러면 Xr의 값이 f(Xr)와 동일한 부호를 갖는 함수 값을 산출하는 2개의 초기 추측 중 하나를 대체한다. 

추가적으로 Bisection 과 False Position 방법을 비교하였을 때, Bisection은 함수의 모양을 고려하지 않는다는 것을 알 수 있다. 이는 함수에 따라 좋을 수도 있고 나쁠 수도 있다.