matlab求有根区间,用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

%     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 著