aitken.m

matlab算法集 matlab算法集

function [x,k] = aitken (x0,b,tol,m,ex,f,dm)
%----------------------------------------------------------------
% Usage:       [x,k] = aitken (x0,b,tol,m,ex,f,dm)
%
% Description: Apply the contraction mapping method with or 
%              without Aitken extrapolation to find a root of:
%
%                  f(x) = 0
%
% Inputs:      x0  = initial guess 
%              b   = gain which ensures that g(x) = x - bf(x) is
%                    a contration.  The gain b must satisfy:
%
%                        0 < bf'(x) < 2
%
%              tol = error tolerance used to terminate search (tol >= 0)
%              m   = maximum number of iterations (m >= 1)
%              ex  = extrapolation switch:
%
%                       0 = contraction mapping
%                       1 = contraction mapping with Aitken extrapolation
%
%              f   = string containing name of user-supplied function 
%                    whose root is to found.  The form of f is:
%
%                       function y = f(x);
%
%                    When f is called, it must return the value f(x).
%              dm  = option display mode. If present, 
%                    intermediate results are displayed.
%
% Outputs:     x = estimated root of f(x)
%              k = number of iterions.  If k < m, then the following
%                  converence criterion was satisfied:
%
%                  |f(x)| < tol 
%----------------------------------------------------------------

% Initialize

   tol = args (tol,0,tol,3,'aitken');
   m   = args (m,1,m,4,'aitken');
   ex  = args (ex,0,1,5,'aitken');
   display = nargin > 6;
   k = 0; 
   x1 = 0;
   x2 = x0;
   x3 = 0;
   f2 = feval(f,x2);
   if display
      fprintf ('\n 0 & %10.7f & %10.7f \\\\',x0,abs(feval(f,x0)));
      fprintf ('\n \\hline');
   end
      
% Find root 

   y = tol + 1;
   while (y > tol) & (k <= m) 
      x3 = x2 - b*f2;
      x0 = x1;
      x1 = x2;
      x2 = x3;
      f2 = feval(f,x3);

      if ex & (k > 2) 
       	 x4 = (x0*x2 - x1*x1)/(x2 - 2*x1 + x0);
         y = abs(feval(f,x4));
         if display
            ffprintf ('\n%2g & %10.7f & %10.7f \\\\',k,x4,y); 
         end
      else 
         y = abs(f2);
         if display
            fprintf ('\n%2g & %10.7f & %10.7f \\\\',k,x2,y); 
         end
      end
      k = k + 1;
   end

% Finalize 

   if ex
      x = x4;
   else
      x = x2;
   end