newton.m

matlab算法集 matlab算法集

function [x,k] = newton (x0,tol,m,f,dm)
%----------------------------------------------------------------------
% Usage:       [x,k] = newton (x0,tol,m,f,dm)
%
% Description: Apply Newton's method to find a root of the n dimensional
%              nonlinear algebraic system: 
%
%                 f(x) = 0
%   
% Inputs:      x0  = n by 1 vector containing initial guess
%              tol = error tolerance used to terminate search (eps >= 0)
%              m   = maximum number of iterations (m >= 1)
%              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 takes the n by 1 vector x
%                    and computes the n by 1 vector y = f(x). Be
%                    sure that y is a column vector, not a row.
%              dm  = optional display mode.  If present, display
%                    intermediate results 
%
% Outputs:     x = n by 1 vector containing estimated root of f(x). 
%              k = the number of iterations performed.  If it is less 
%                  than the user-specified maximum, m, then the 
%                  following converence criterion was satisfied:
%
%                      ||f(x)|| < tol 
%----------------------------------------------------------------------

% Initialize 

   chkvec (x0,1,'newton');
   tol = args (tol,0,tol,2,'newton');
   m = args (m,1,m,3,'newton');
   chkfun (feval(f,x0),4,'newton')

   display = nargin > 4;
   k = 1;
   x = x0;
   n = length (x);
   y  = zeros (n,1);
   dx = zeros (n,1);
   J  = zeros (n,n);
   y = feval(f,x);

   if display
      fprintf ('\n%2g & ',0);
      for i = 1 : n
         fprintf ('%10.7f & ',x(i));
      end 
      fprintf ('%10.7f \\\\',norm(y,inf)); 
      fprintf ('\n \\hline');
   end
   
% Find root 

  r = 2*tol;
  while (r > tol) & (k <= m)
     J = jacobian (x,f);
     dx = gauss (J,y);
     x = x - dx;
     y = feval(f,x);
     r = norm (y,inf);
     if display
        printf ('\n%2g & ',k);
        for i = 1 : n
           printf ('%10.7f & ',x(i));
        end  
        printf ('%10.7f \\\\',r);
     end         
     k = k + 1;
  end
  k = k - 1;
%----------------------------------------------------------------------