muller.m

matlab算法集 matlab算法集

function [x,k] = muller (x0,x1,x2,tol,m,f,dm)
%----------------------------------------------------------------
% Usage:       [x,k] = muller (x0,x1,x2,tol,m,f)
%
% Description: Apply the Muller method to find a root of:
%
%                   f(x) = 0
%
% Inputs:       x0  = first initial guess
%               x1  = second initial guess (x1 <> x0)
%               x2  = third initial guess (x2 <> x0, x2 <> x1)
%               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
%----------------------------------------------------------------

% Initialize 

   if (abs(x1-x0) < eps) | (abs(x2-x0) < eps) | (abs(x2-x1) < eps)
      fprintf ('Secant requires that x2 <> x1 <> x0.\n');
      return
   end
   tol = args (tol,0,tol,4,'muller');
   m   = args (m,1,m,5,'muller');
   display = nargin > 6;

   f0 = feval(f,x0);
   f1 = feval(f,x1);
   f2 = feval(f,x2);
   if display
      fprintf ('\n 0 & %10.7f %+10.7fj & %10.7f \\\\',x0.x,x0.y,cmag(f0));
      fprintf ('\n 1 & %10.7f %+10.7fj & %10.7f \\\\',x1.x,x1.y,cmag(f1));
      fprintf ('\n 2 & %10.7f %+10.7fj & %10.7f \\\\',x2.x,x2.y,cmag(f2));
      fprintf ('\n \\hline');
   end

% Find root 
   
   k = 1;
   y = eps + 1;
   while (y > eps) & (k <= m)

      h0 = x0 - x2;
      h1 = x1 - x2;
      d = h0*h1*(h1-h0);
      c0 = f2;

      ca = h1*h1*(f0-f2);
      cb = h0*h0*(f1-f2);
      c1 = (ca-cb)/d;

      ca = h0*(f1-f2);
      cb = h1*(f0-f2);
      c2 = (ca - cb)/d;

      c = 4*c0;
      c = sqrt(c1*c1 - c*c2);
      ca = c1 + c;
      cb = c1 - c;
      if abs(cb) > abs(ca) 
         ca = cb;
      end
      c = 2*c0;
      cb = c/ca;
      x3 = x2 - cb;
      
      x0 = x1;
      f0 = f1;
      x1 = x2;
      f1 = f2;
      x2 = x3;
      f2 = feval(f,x3);
      y = abs(f2);
      if display
         fprintf ('\n%2g & %10.7f %+10.7fj & %10.7f \\\\',k+2,x2.x,x2.y,y);
      end
      k = k + 1;
   end

% Finalize 

   x = x2;   
   k = k-1;