Home > K25104 > NonlinearSystems > secant_method.m

secant_method

PURPOSE ^

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

SYNOPSIS ^

function [ x,k ] = secant_method( f,a,b,tol,max )

DESCRIPTION ^

 Implements secant method algorithm to find root of f(x)=0
 Input arguments:
   f, function handle or vector of polynomial coeffiencts
   a,b, end points of initial search interval
   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 ] = secant_method( f,a,b,tol,max )
0002 % Implements secant method algorithm to find root of f(x)=0
0003 % Input arguments:
0004 %   f, function handle or vector of polynomial coeffiencts
0005 %   a,b, end points of initial search interval
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 a>=b;    
0014     error('specifed interval invalid');
0015 elseif tol<=0;
0016     error('tol must be positive');
0017 elseif max<=0;
0018     error('max must be positive');
0019 elseif isa(f,'function_handle');
0020     % do nothing
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 function handle 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         f=matlabFunction(v);    % this converts the expression to a
0030                                 % function handle
0031     end
0032 else
0033     error(['input must be a function handle or vector holding ',...
0034             'polynomial coefficients']);
0035 end
0036 
0037 % start of algorithm
0038 fa=feval(f,a);  % evaluate f(a)
0039 fb=feval(f,b);  % evaluate f(b)
0040 k=0;            % initialize counter
0041 
0042 if fa==0;           % then a is a root
0043     x=a; return;
0044 elseif fb==0;       % then b is a root
0045     x=b; return;    
0046 end
0047 fc=fb;              % initialise f_c to start loop
0048 
0049 while abs(fc)>tol && k<=max;    
0050     
0051     denom=fb-fa;
0052     if abs(denom)<1e-12; % check size of denominator
0053         disp(['denominator gets very close to zero, possible'...
0054             'loss of significance'])    % display warning
0055     end
0056     
0057     c=b-fb*((b-a)/(denom));     % calculate next iterate
0058     
0059     fc=feval(f,c);              % evaluate f(c)
0060     a=b;        % update a
0061     b=c;        % update b
0062     fa=fb;      % update f(a)
0063     fb=fc;      % update f(b)
0064    
0065     k=k+1;      % increment counter
0066 end
0067 
0068 x=c;    % set x to the most recent iterate
0069 
0070 if k==max+1;
0071     disp('maximum number of iterations reached');
0072 end
0073 
0074 end

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