findzero.m

寻找0点或者方程求解寻找0点或者方程求解寻找0点或者方程求解

function [ksi,fval,iter] = findzero(hFcn,x0,TOL,varargin)
%findzero  Scalar nonlinear zero finding. 
%   x = FINDZERO(hFcn,x0) tries to find a zero of the function hFcn near x0, 
%   if x0 is a scalar.  It first finds an interval containing x0 where the 
%   function values of the interval endpoints differ in sign, then searches 
%   that interval for a zero. hFcn accepts real scalar input x and returns 
%   a real scalar function value F, evaluated at x.  The value x returned 
%   by FINDZERO is near a point where hFcn changes sign (if hFcn is continuous), 
%   or NaN if the search fails.  
%
%   x = FINDZERO(hFcn,x0), where x0 is a vector of length 2, assumes x0 is an 
%   interval where the sign of hFcn(x0(1)) differs from the sign of hFcn(x0(2)).
%   An error occurs if this is not true.  Calling FINDZERO with an interval  
%   guarantees FINDZERO will return a value near a point where hFcn changes 
%   sign.
%
%   x = FINDZERO(hFcn,x0), where x0 is a scalar value, uses x0 as a starting 
%   guess. FINDZERO looks for an interval containing a sign change for hFcn and 
%   containing x0.  If no such interval is found, NaN is returned.  
%   In this case, the search terminates when the search interval 
%   is expanded until an Inf, NaN, or complex value is found. 
%
%   Examples
%     hFcn can be specified using @:
%        x = fzero(@sin,3)
%     returns pi.
%        x = fzero(@sin,3,optimset('disp','iter')) 
%     returns pi, uses the default tolerance and displays iteration information.
%
%     hFcn can also be an anonymous function:
%        x = fzero(@(x) sin(3*x),2)
%
%   If hFcn is parameterized, you can use anonymous functions to capture the 
%   problem-dependent parameters. Suppose you want to solve the equation given
%   in the function MYFUN, which is parameterized by its second argument A. Here
%   MYFUN is an M-file function such as
%
%       function f = myfun(x,a)
%       f = cos(a*x);
%
%   To solve the equation for a specific value of A, first assign the value to A.
%   Then create a one-argument anonymous function that captures that value of A 
%   and calls MYFUN with two arguments. Finally, pass this anonymous function to 
%   FINDZERO:
%
%       a = 2; % define parameter first
%       x = fzero(@(x) myfun(x,a),0.1)
%   
%   Limitations
%        x = fzero(@(x) abs(x)+1, 1) 
%     returns NaN since this function does not change sign anywhere on the 
%     real axis (and does not have a zero as well).
%        x = fzero(@tan,2)
%     returns x near 1.5708 because the discontinuity of this function near the 
%     point x gives the appearance (numerically) that the function changes sign at x.
%

%  This algorithm was originated by T. Dekker.  An Algol 60 version,
%  with some improvements, is given by Richard Brent in "Algorithms for
%  Minimization Without Derivatives", Prentice-Hall, 1973.  A Fortran
%  version is in Forsythe, Malcolm and Moler, "Computer Methods
%  for Mathematical Computations", Prentice-Hall, 1976.

% 初始化
iter = 0;
ksi = NaN;
fval = NaN;
if nargin < 3, TOL = eps; end
% 输入值x0是区间 x0 = [a,b]
if (length(x0) == 2) 
    xa = x0(1); xa0=xa;
    xb = x0(2); xb0=xb;
    ya = feval(hFcn,xa,varargin{:});
    yb = feval(hFcn,xb,varargin{:});
    if any(~isfinite([ya yb])) || any(~isreal([ya yb]))
        disp('区间端点的函数值必须是有限的实数.')
        return
    end
    ya0 = ya; yb0 = yb;
    
    if ( ya == 0 )
        ksi = xa; fval = ya;
        return
    elseif ( yb == 0)
        ksi = xb; fval = yb;
        return
    elseif (ya > 0) == (yb > 0)
        disp('区间端点的函数值必须不同号.')
        return
    end    
% 输入值x0是初始试探零点
elseif (length(x0) == 1)
    fx = feval(hFcn,x0,varargin{:});
    if fx == 0
        ksi = x0; fval = fx;
        return
    elseif ~isfinite(fx) || ~isreal(fx)
        disp('区间端点的函数值必须是有限的实数.');
        return
    end
    if x0 ~= 0, 
        dx = x0/50;
    else 
        dx = 1/50;
    end
    % 找出函数值正负号改变的位置.
    twosqrt = sqrt(2); 
    xa = x0; ya = fx; xb = x0; yb = fx;
    while (ya > 0) == (yb > 0)
        dx = twosqrt*dx;
        xa = x0 - dx;  
        ya = feval(hFcn,xa,varargin{:});
        if ~isfinite(ya) || ~isreal(ya)
            ksi = NaN; fval = NaN;
            return
        end
        if (ya > 0) ~= (yb > 0) 
            break   % 找到了函数值变号的区间
        end
        xb = x0 + dx;  yb = feval(hFcn,xb,varargin{:});
        if ~isfinite(yb) || ~isreal(yb)
            ksi = NaN; fval = NaN;
            return
        end
    end % while
    xa0 = xa; ya0 = ya; xb0 = xb; yb0 = yb;
else
    disp('第二个输入宗量必须是长度为1或2的数组.');
    return
end % if (length(x0) == 2

yc = yb;
% 主循环. 一旦中途找到根，程序就退出.
while yb ~= 0
    % 保证xb是迄今最好的近似, xa是xb以前的值, 在xc与xb处f(x)正负号相反
    if (yb > 0) == (yc > 0)
        xc = xa;  yc = ya;
        dx = xb - xa;  er = dx;
    end
    if abs(yc) < abs(yb)
        xa = xb;  xb = xc;  xc = xa;
        ya = yb;  yb = yc;  yc = ya;
    end
    % 收敛性检测
    del = 0.5*(xc - xb);
    xtol = 2.0*TOL*max(abs(xb),1.0);
    if (abs(del) <= xtol) || (yb == 0.0) 
        break 
    end
    % 二分法或插值法的选择
    if (abs(er) < xtol) || (abs(ya) <= abs(yb))
        dx = del;  er = del;    % 选择二分法
    else
        % 以下是插值方法
        s = yb/ya;
        if (xa == xc)   % 线性插值
            p = 2.0*del*s;
            q = 1.0 - s;
        else    % 反二次插值
            q = ya/yc;
            r = yb/yc;
            p = s*(2.0*del*q*(q - r) - (xb - xa)*(r - 1.0));
            q = (q - 1.0)*(r - 1.0)*(s - 1.0);
        end;
        if p > 0, q = -q; else p = -p; end;
        % 判别是否接受插值点
        if (2.0*p < 3.0*del*q - abs(xtol*q)) && (p < abs(0.5*er*q))
            er = dx;  dx = p/q;
        else
            dx = del;  er = del;
        end;
    end % 结束插值
    % 预备下一个迭代
    xa = xb;
    ya = yb;
    if abs(dx) > xtol, xb = xb + dx;
    elseif xb > xc, xb = xb - xtol;
    else xb = xb + xtol;
    end
    yb = feval(hFcn,xb,varargin{:});
    iter = iter + 1;
end % 结束主循环
ksi = xb;
fval = feval(hFcn,ksi,varargin{:});
if abs(fval) <= max(abs(ya0),abs(yb0))
    msg = sprintf('在区间 [%g, %g] 内找到零点',xa0,xb0);
else
    msg = sprintf('区间 [%g, %g] 内可能包含f(x)的奇异点',xa0,xb0);
end
disp(msg)