用matlab得到有根区间向量

function Xb=brackPlot(fun,xmin,xmax,nx)
% BrackPlot find subintervals on x that contain sign changes of f(x)
% Synopsis: Xb=brackPlot(fun,xmin,xmax,nx)
%           Xb=brackPlot(fun,xmin,xmax)
% Input:    fun=(string) name of mfile function that evaluates f(x)
%           xmin,xmax=endpoints of interval to subdivide into brackets
%           nx=(optional) number of samples along x axis used to test for
%           brackets.The interval xmin<=x<=xmax is divided into nx-1
%           subinterval. Default:nx=20
%  Output:  Xb=two column matrix of bracket limits.
%           Xb(k,1) is the left (lower x value) bracket and
%           Xb(k,2) is the ringht bracket for the k^th potenntial root.
%           If no bracket are found,Xb=[].
if nargin < 4
    nx=20;
end

% ----Plot f(x) on interval xmin <= x<= xmax
xp=linspace(xmin,xmax);
yp=feval(fun,xp);
plot(xp,yp,[floor(xmin) ceil(xmax)],[0,0]);
grid on;
xlabel('x');
ylabel(['f(x) defined in ',fun,'.m']);

%----save data used to draw boxes that indicate brackets
ytop=max(yp);ybot=min(yp);              % y coordinates of the box
ybox=0.05*[ybot ytop ytop ybot ybot];   % around a bracket
c=[0.7 0.7 0.7];                        % RGB color used to fill the box

%----begin search for brackets
x=linspace(xmin,xmax,nx);              % vector of potential bracket limits
f=feval(fun,x);                        % vector of f(x) values at potential brackets
nb=0;  Xb=[];                          % Xb is null unless brackets are found
for k=1:length(f)-1
    if sign(f(k))~=sign(f(k+1))        % True if sign of f(x) changes in the interval
        nb=nb+1;
        Xb(nb,1)=x(k);
        Xb(nb,2)=x(k+1);              % save left and right ends of the brackets
        hold on;
        fill([x(k) x(k) x(k+1) x(k+1) x(k)],ybox,c);
                                      % add filled box
    end
end
hold off
if isempty(Xb)                        % free advice
    warning('No brackets found.Check [xmin,xmax] or increase nx');
end
                                 
%用matlab得到有根区间向量以便其他求根程序使用
%create date:2008.3.25
%  区间划分法算法：
%  given:f(x),xmin,xmax,n
%  dx=(xmax-xmin)/n
%  a=xmin
%  i=0
%  while i<n
%     i=i+1
%     b=a+dx
%     if f(x) change sign in [a,b]
%         save [a,b] for further root-finding
%     end
%     a=b
%  end

%  在测试部分“if f(x) change sign in [a,b]”，简单的实现方法为：
%       f(a)*f(b)< 0
%  这个方法在浮点数运算时会存在问题，如果f(a) f(b)的值都很小（当区间的两端都很接近某个根时），
%  上面两个数的乘积可能会导致下溢从而导致对符号变化的检测出错。如：
%   format long e
%   fa=1e-120;  fb=-2e-300
%   fa*fb=0
%  所以，本函数中采用的是sign函数进行判别
%  不管fa,fb值多小，只要它们符号不同，表达式 sign(fa)~=sign(fb)总会返回真实值
%  资料来源：数值方法和matlab实现与应用  机械工业出版社  (美)Gerald Recktenwald 著