Home > K25104 > NonlinearSystems > Newton_Raphson.m

Newton_Raphson

PURPOSE ^

Implements Newton method to find root of f(x)=0

SYNOPSIS ^

function [ x,k ] = Newton_Raphson( f,initial,tol,max )

DESCRIPTION ^

 Implements Newton method to find root of f(x)=0
 Input arguments:
   f, function handle or vector of polynomial coefficients
   initial, initial iterate
   tol, tolerance
   max, maximum number of iterations
 Output arguments:
   x, solution
   k, number of iterations to convergence

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ x,k ] = Newton_Raphson( f,initial,tol,max )
0002 % Implements Newton method to find root of f(x)=0
0003 % Input arguments:
0004 %   f, function handle or vector of polynomial coefficients
0005 %   initial, initial iterate
0006 %   tol, tolerance
0007 %   max, maximum number of iterations
0008 % Output arguments:
0009 %   x, solution
0010 %   k, number of iterations to convergence
0011 
0012 % first check user inputs
0013 if tol<=0;
0014     error('tolerance must be positive');
0015 elseif max<=0;
0016     error('maximum number of iterations must be positive');
0017 elseif isa(f,'sym');
0018     df = diff(f);
0019     f = matlabFunction(f);
0020     df = matlabFunction(df);
0021 elseif isa(f,'double'); 
0022     [n,m]=size(f);  % find the size of f
0023     if n~=1 && m~=1; % then f is not a vector
0024         error(['input must be a symbolic function or vector holding ',...
0025             'polynomial coefficients']);
0026     else 
0027         v=poly2sym(f);  % this converts the vector of coefficients to an
0028                         % expression for the polynomial
0029         df = diff(v);
0030         f=matlabFunction(v);    % this converts the expression to a
0031                                 % function handle
0032         df = matlabFunction(df);
0033     end
0034 else
0035     error(['input must be a symbolic function or vector holding ',...
0036             'polynomial coefficients']);
0037 end
0038 
0039 % start of algorithm
0040 x=zeros(max,1);         % initialise vector to hold iterates
0041 x(1)=initial;           % initialize x(1)
0042 fx=feval(f,x(1));       % evaluate f(x)
0043 divx=feval(df,x(1));   % evaluate f'(x)
0044 k=2;                % initialize iteration counter
0045 l=0;    % initialize another counter to check if iterates are diverging
0046 d0=0;   % initialize d0 to zero
0047 
0048 while abs(fx)>tol && k<=max   
0049     
0050     x(k)=x(k-1)-fx/divx;    % update our estimate of the root
0051     
0052     test=find(x(1:k-1)==x(k));
0053     if isempty(test)==0;    % if our estimate of the root matches any of
0054                             % the previous estimates then we fall into
0055                             % an infinite loop
0056         disp('algorithm fell into infinite loop');
0057         x=x(1:k); return;   % return vector of iterates so far
0058     end
0059     
0060     d1=abs(x(k)-x(k-1));    % calculate d1
0061     if d1>d0;   % then distance between succesive iterates is increasing
0062         l=l+1;  % increment l
0063         if l>20;
0064             disp('distance between succesive iterates is increasing');
0065             x=x(1:k); return;   % return vector of iterates so far
0066         end
0067     else
0068         l=0;    % reset l
0069     end
0070     d0=d1;      % update d0
0071 
0072     fx=feval(f,x(k));       % evaluate f(x)
0073     divx=feval(df,x(k));   % evaluate f'(x)
0074     
0075     if abs(divx)<1e-12;
0076         disp('denominator vanishes'); 
0077         x=x(1:k); return;   % return vector of iterates so far
0078     end
0079     
0080     k=k+1;  % increment counter
0081 end
0082 
0083 if k==max+1;
0084     disp('maximum number of iterations reached');
0085 end
0086 x=x(k-1);    % set x to the most recent iterate
0087 
0088 end

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005