function [root,success] =bisect_method (a, b, eps, f)
% Accompanying program for the text
%
% Classical and Modern Numerical Analysis:
% Theory, Methods and Practice
% by Azmy S. Ackleh, Edward J. Allen,
% R. Baker Kearfott, and Padmanabhan Seshaiyer
%
% [root, success] = bisect_method (a, b, eps, func) returns the result
% of the method of bisection, with starting interval [a, b],
% tolerance eps, and with function defined by y = f(x). For example,
% suppose an m-file xsqm2.m is available in Matlab's working
% directory, with the following contents:
% function [y] = xsqm2(x)
% y = x^2-2;
% return
% Then, issuing
% [root,success] = bisect_method (1, 2, 1e-10, 'xsqm2')
% from Matlab's command window will cause an approximation to
% the square root of 2 that, in the absence of excessive roundoff
% error, has absolute error of at most 10^{-16}
% to be returned in the variable root, and success to be set to
% 'true'.
%
% success is set to 'false' if f(a) and f(b) do not have the same
% sign. success is also set to 'false' if the tolerance cannot be met.
% In either case, a message is printed, and the midpoint of the present
% interval is returned in the variable root.
% Algorithm 2.1 is basically used. However, Theorem 2.2 is used to
% compute the number of iterations N required for bisection to achieve
% that tolerance. This is to avoid an infinite loop that might be
% caused by an excessively small tolerance (in which case the midpoint of
% the interval is repeatedly rounded to one of the end points and the
% tolerance is never met).
error=b-a;
fa=feval(f,a);
fb=feval(f,b);
success = true;
% First, handle incorrect arguments --
if (fa*fb > 0)
disp('Error: f(a)*f(b)>0');
success = false;
root = a + (b-a)/2;
return
end
if (eps <=0)
disp('Error: eps is less than or equal to 0')
success = false;
root = a + (b-a)/2;
return
end
if (b < a)
disp('Error: b < a')
success = false;
root = (a+b)/2;
return
end
% Set N to be the smallest integer such that N iterations of bisection
% suffices to meet the tolerance --
N = ceil( log((b-a)/eps)/log(2) - 1 )
% This is where we actually do Algorithm 2.1 --
disp(' -----------------------------');
disp(' Error Estimate ');
disp(' -----------------------------');
for i=1:N
x= a + (b-a)/2;
fx=feval(f,x);
if(fx*fa > 0)
a=x;
else
b=x;
end
error=b-a;
disp(sprintf(' %12.4e %12.4e', error, x));
end
% Finally, check to see if the tolerance was actually met. (With
% additional analysis of the minimum possible relative error as in Theorem
% 1.6 and Definition 1.5, unreasonable values of epsilon for the particular
% floating point system can be determined before the loop on i, avoiding
% unnecessary work.)
error = (b-a)/2;
root = a + (b-a)/2;
if (error > eps)
disp('Error: epsilon is too small for tolerance to be met');
success = false;
return
end