bisect.m

matlab算法集 matlab算法集

function [x,k] = bisect (x0,x1,tol,m,f,dm)
%----------------------------------------------------------------
% Usage:       [x,k] = bisect (x0,x1,tol,m,f,dm)
%
% Description: Apply the bisection method to find a root of:
%
%                   f(x) = 0
%
% Inputs:       x0  = lower limit of search
%               x1  = upper limit of search (x1 > x0)
%               tol = error tolerance used to terminate search 
%                     (tol >= 0)
%               m   = maximum number of iterations (m >= 1)
%               f   = string containing name of user-supplied
%                     function.  The function f is of the form:
%
%                     function y = f(x)
%
%               dm  = optional display mode.  If present,
%                     intermediate results are displayed.
%
% Outputs:      x = Estimated root of f(x) = 0 
%               k = number of iterations peroformed. If k < m,
%                   then the following convergence criterion
%                   was satisfied:
%
%                   |f(x)| < eps
%   
% Notes:        1. x0 and x1 must be selected such that 
%                  f(x0)*f(x1) < 0 to ensure that there is a 
%                  root in [x0,x1]
%
%               2. bisect should always converge for a continuous
%                  f if tol and m are sufficiently large.  
%----------------------------------------------------------------

% Initialize 

   x1  = args (x1,x0,x1,2,'bisect');
   tol = args (tol,0,tol,3,'bisect');
   m   = args (m,1,m,4,'bisect');
   display = nargin > 5;
   f0 = feval (f,x0);
   f2 = feval (f,x1);
   if f0*f2 >= 0
      fprintf ('In bisect it is required that f(x0)f(x1) < 0.\n');
      return
   end

% Find root 

   k = 1;
   y = tol + 1;
   while (y > tol) & (k <= m)
      x2 = (x0 + x1)/2;
      f2 = feval(f,x2);
      if f0*f2 < 0
         x1 = x2;
      else 
         x0 = x2;
         f0 = f2;
      end
      y = abs(f2);
      if display
         fprintf ('\n%2i & %10.7f & %10.7f & %10.7f \\\\',k,x2,x1-x0,y);
      end
      k = k + 1;	   
   end

% Finalize 

   x = x2;
   k = k - 1;