% Matlab code bisect.m
% For "Applied Numerical Linear Algebra", Algorithm 1.1 and Question 1.21
% James Demmel, Aug 1995, modified May 1997
%
% function x = bisect(coeff,a,b,tol)
%
% Find a zero of a polynomial p(x)=0
% given an interval [a,b] where p(a)*p(b) < 0.
%
% Inputs:
% coeff = vector of polynomial coefficients
% p(x) = sum from i=0 to n of coeff(i+1)*x^(n+1-i)
% where n = length(coeff)+1
% a,b = interval containing (at least) one zero
% tol = stopping tolerance. stop bisecting if
% the interval containing a zero has width < 2*tol
%
% Outputs:
% x = [xlow,xhigh] = interval containing a zero
%
function x = bisect(coeff,a,b,tol)
xlow = a;
xhigh = b;
plow = polyval(coeff,xlow);
phigh = polyval(coeff,xhigh);
while xhigh - xlow > 2*tol,
xmid = (xlow + xhigh)/2;
pmid = polyval(coeff,xmid);
if pmid*plow < 0,
xhigh = xmid;
phigh = pmid;
elseif pmid*phigh < 0,
xlow = xmid;
plow = pmid;
else
xlow = xmid;
xhigh = xmid;
end
end
x = [xlow, xhigh];