Bisection 방법으로 근을 구하는 Matlab 함수

기계시스템디자인공학과 김성환 교수

(1) [eqn: f(m)= sqrt {{9.81 TIMES m} over {0.25}} tanh LEFT ( sqrt {{9.81 TIMES 0.25} over {m}} TIMES 4 RIGHT ) -36] 의 근 구하기 위해

근을 찾을 구간을 정한다. 하한경계값을 [eqn: x _{l} =40], 상한경계값을 [eqn: x _{u} =200]이라 하면

[eqn: f(40)=] -5.9870, [eqn: f(200)=] 0.8603 이므로 (40, 200) 구간에 해가 존재한다.

>> fm = @(m) sqrt(9.81*m/0.25)*tanh(sqrt(9.81*0.25/m)*4)-36

>> fm(40)

>> fm(200)

(2) 그래프를 그려 개형을 알아본다.

- 근은 대략 145 부근에 있음

>> fplot(fm, [40,200]); grid

(3) 알고리즘

- [eqn: x _{r} =(x _{l} +x _{u} )/2] 를 위한다.

[eqn: f(x _{l} ) BULLET f(x _{r} )<0`] 이면 [eqn: x _{r}]을 새로운 [eqn: x _{u}]로 삼음

[eqn: f(x _{r} ) BULLET f(x _{u} )<0`] 이면 [eqn: x _{r}]을 새로운 [eqn: x _{l}]로 삼음

상대근사오차 [eqn: epsilon _{a} = LEFT | {x _{r} ^{"new"} -x _{r} ^{old}} over {x _{r} ^{"new"}} RIGHT | TIMES 100] 이 허용오차 [eqn: epsilon _{s}]보다 작아지면 반복 종료

(4) 함수값 계산을 코드에 포함

function root = bsect(xl, xu)

xr = xl;

while(1)

xrold = xr;

xr = (xl + xu) / 2;

fxl = sqrt(9.81*xl/0.25)*tanh(sqrt(9.81*0.25/xl)*4)-36;

fxr = sqrt(9.81*xr/0.25)*tanh(sqrt(9.81*0.25/xr)*4)-36;

test = fxl * fxr;

if test < 0

xu = xr;

elseif test > 0

xl = xr;

else

ea = 0;

end

ea = abs( (xr-xrold)/xr ) * 100;

if ea <= 0.0001; break; end

end

root = xr;

end

>> bsect(40, 200)

(5) 반복 횟수를 리턴

function [root, niter] = bsect(xl, xu)

niter = 0;

xr = xl;

while(1)

xrold = xr;

xr = (xl + xu) / 2;

fxl = sqrt(9.81*xl/0.25)*tanh(sqrt(9.81*0.25/xl)*4)-36;

fxr = sqrt(9.81*xr/0.25)*tanh(sqrt(9.81*0.25/xr)*4)-36;

test = fxl * fxr;

if test < 0

xu = xr;

elseif test > 0

xl = xr;

else

ea = 0;

end

niter = niter + 1;

ea = abs( (xr-xrold)/xr ) * 100;

if ea <= 0.0001; break; end

end

root = xr;

end

>> [root, niter] = bsect(40, 200)

(6) 함수값 계산을 부함수로 - 내포화

function [root, niter] = bsect(xl, xu)

niter = 0;

xr = xl;

while(1)

xrold = xr;

xr = (xl + xu) / 2;

test = func(xl) * func(xr);

if test < 0

xu = xr;

elseif test > 0

xl = xr;

else

ea = 0;

end

niter = niter + 1;

ea = abs( (xr-xrold)/xr ) * 100;

if ea <= 0.0001; break; end

end

root = xr;

end

function f = func(x)

f = sqrt(9.81*x/0.25)*tanh(sqrt(9.81*0.25/x)*4)-36;

end

>> [root, niter] = bsect(40, 200)

(7) 함수를 외부에서 인수로 입력

function [root, niter] = bsect(func, xl, xu)

niter = 0;

xr = xl;

while(1)

xrold = xr;

xr = (xl + xu) / 2;

test = func(xl) * func(xr);

if test < 0

xu = xr;

elseif test > 0

xl = xr;

else

ea = 0;

end

niter = niter + 1;

ea = abs( (xr-xrold)/xr ) * 100;

if ea <= 0.0001; break; end

end

root = xr;

end

>> func = @(m) sqrt(9.81*m/0.25)*tanh(sqrt(9.81*0.25/m)*4)-36

>> [root, niter] = bsect(func, 40, 200)

▢