Home > K25104 > NonlinearSystems > inv_quad_it.m

inv_quad_it

PURPOSE ^

Implements inverse quadratic interpolation method to find root of f(x)=0

SYNOPSIS ^

function [ x,k ] = inv_quad_it( f,x0,x1,x2,tol,max )

DESCRIPTION ^

 Implements inverse quadratic interpolation method to find root of f(x)=0
 Input arguments:
   f, function handle or vector of polynomial coefficients
   x0,x1,x2, initial iterates
   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 ] = inv_quad_it( f,x0,x1,x2,tol,max )
0002 % Implements inverse quadratic interpolation method to find root of f(x)=0
0003 % Input arguments:
0004 %   f, function handle or vector of polynomial coefficients
0005 %   x0,x1,x2, initial iterates
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,'function_handle');
0018     % do nothing
0019 elseif isa(f,'double'); 
0020     [n,m]=size(f);  % find the size of f
0021     if n~=1 && m~=1; % then f is not a vector
0022         error(['input must be a function handle or vector holding ',...
0023             'polynomial coefficients']);
0024     else 
0025         v=poly2sym(f);  % this converts the vector of coefficients to an
0026                         % expression for the polynomial
0027         f=matlabFunction(v);    % this converts the expression to a
0028                                 % function handle
0029     end
0030 else
0031     error(['f must be a function handle or vector holding ',...
0032             'polynomial coefficients']);
0033 end
0034 
0035 % start of algorithm
0036 x=zeros(max,1);     % initialize vector to hold iterates
0037 x(1)=x0;x(2)=x1;x(3)=x2;      % initialize first three iterates
0038 f0=feval(f,x(1));           % evaluate f(x0)
0039 f1=feval(f,x(2));           % evaluate f(x1)
0040 f2=feval(f,x(3));           % evaluate f(x2)
0041 k=4;                % initialize iteration counter
0042 l=0;    % initialize another counter to check if iterates are diverging
0043 d0=abs(x2-x1);      % initialize d0 to hold the distance between the
0044                     % last two iterates
0045 
0046 while abs(f2)>tol && k<=max;   
0047     
0048     if abs(f0-f1)<1e-12;
0049         disp('denominator vanishes'); 
0050         x=x(1:k); return;     % return vector of iterates so far
0051     elseif abs(f1-f2)<1e-12;
0052         disp('denominator vanishes'); 
0053         x=x(1:k); return;     % return vector of iterates so far
0054     elseif abs(f0-f2)<1e-12;
0055         disp('denominator vanishes'); 
0056         x=x(1:k); return;     % return vector of iterates so far
0057     end
0058     
0059     % update our estimate of the root
0060     x(k)=x0*f1*f2/((f0-f1)*(f0-f2))+f0*x1*f2/((f1-f0)*(f1-f2))+...
0061         f0*f1*x2/((f2-f0)*(f2-f1));   
0062     
0063     d1=abs(x(k)-x(k-1));    % update d1
0064     if d1>d0;   % then distance between succesive iterates is increasing
0065         l=l+1;  % increment l
0066         if l>max/10;    % check if the distance between successive iterates
0067                         % has been increasing for the last max/10
0068                         % iterations
0069             disp('distance between succesive iterates is increasing');
0070             x=x(1:k); return;   % return vector of iterates so far
0071         end
0072     else
0073         l=0;    % reset l
0074     end
0075     d0=d1;      % update d0
0076 
0077     x0 = x1;
0078     x1 = x2;
0079     x2 = x(k);
0080     f0=f1;              % update f0
0081     f1=f2;              % update f1
0082     f2=feval(f,x(k)); % update f2
0083     
0084     k=k+1;  % increment counter
0085 end
0086 
0087 if k==max+1;
0088     disp('maximum number of iterations reached');
0089 end
0090 x=x(k-1);    % set x to the most recent iterate
0091 
0092 end

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