Home > K25104 > NonlinearSystems > Halley.m

Halley

PURPOSE ^

Implements Halley's method to find root of f(x)=0

SYNOPSIS ^

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

DESCRIPTION ^

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

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