Home > K25104 > NonlinearSystems > Bus_Dekker.m

Bus_Dekker

PURPOSE ^

Implements the Bus and Dekker algorithm to find a solution of f(x)=0;

SYNOPSIS ^

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

DESCRIPTION ^

 Implements the Bus and Dekker algorithm to find a solution of f(x)=0;
 Input arguments:
   f, function handle or vector of polynomial coefficients
   b, initial iterate
   a, initial contrapoint
   tol, tolerance
   max, maximum number of iterations
 Output arguments:
   x, solution
   k, number of iterations

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ x,k ] = Bus_Dekker( f,a,b,tol,max )
0002 % Implements the Bus and Dekker algorithm to find a solution of f(x)=0;
0003 % Input arguments:
0004 %   f, function handle or vector of polynomial coefficients
0005 %   b, initial iterate
0006 %   a, initial contrapoint
0007 %   tol, tolerance
0008 %   max, maximum number of iterations
0009 % Output arguments:
0010 %   x, solution
0011 %   k, number of iterations
0012 
0013 % first check user inputs
0014 if a>=b;    
0015     error('specifed interval invalid');
0016 elseif tol<=0;
0017     error('tolerance must be positive');
0018 elseif max<=0;
0019     error('maximum number of iterations must be positive');
0020 elseif isa(f,'function_handle');
0021     % do nothing
0022 elseif isa(f,'double'); 
0023     [n,m]=size(f);  % find the size of f
0024     if n~=1 && m~=1; % then f is not a vector
0025         error(['input must be a function handle or vector holding ',...
0026             'polynomial coefficients']);
0027     else 
0028         v=poly2sym(f);  % this converts the vector of coefficients to an
0029                         % expression for the polynomial
0030         f=matlabFunction(v);    % this converts the expression to a
0031                                 % function handle
0032     end
0033 else
0034     error(['input must be a function handle or vector holding ',...
0035             'polynomial coefficients']);
0036 end
0037 
0038 % start of algorithm
0039 fa=feval(f,a);     % evaluate f(a)
0040 fb=feval(f,b);     % evaluate f(b)
0041 
0042 if fa==0;          % check incase the solution is at one of the endpoints
0043     x=a; return;
0044 elseif fb==0;   
0045     x=b; return;
0046 elseif fa*fb>0;   % if f(a),f(b) have the same sign
0047         error('function values at endpoints have the same sign');
0048 end
0049 
0050 c=a;    % we will need to hold the previous two iterates, so initialize
0051         % another variable. In the first instance set c=a
0052 fc=fa;      % store f(c)
0053 k=0;        % initialize counter
0054 
0055 while abs(b-a)>tol && k<=max;
0056     
0057     m=(a+b)/2;         % calculate the midpoint of the interval
0058     s=b-((b-c)*fb)/(fb-fc);     % calculate the approximation to x
0059                                 % given by secant method
0060     temp=b;    % hold b
0061     
0062     if or((s>b & s<m),(s<b & s>m)); % check if s lies between b and m
0063         b=s;   % update b
0064     else
0065         b=m;   % update b
0066     end
0067     fb=feval(f,b);      % evaluate f(b)
0068     
0069     if fa==0;           % check if a is a root
0070         x=a; return;
0071     elseif fb==0;       % check if b is a root
0072         x=b; return;
0073     elseif fa*fb>0;     %  if f(a),f(b) have the same sign
0074         a=c;            % update a
0075         fa=feval(f,a);  % evaluate f(a)
0076     end
0077     
0078     c=temp;             % update c
0079     fc=feval(f,c);      % evaluate f(c)
0080     
0081     
0082     if abs(fa)<abs(fb); % then we need to interchange a and b, so that
0083                         % the iterate is the best approximation to x
0084         swap=a;         % hold a
0085         a=b; b=swap;    % swap b and a
0086         fa=feval(f,a);  % revaluate f(a)
0087         fb=feval(f,b);  % revaluate f(b)
0088     end
0089    
0090     k=k+1;      % increment k
0091 end
0092 
0093 x=b;    % set x to be the last iterate
0094 
0095 if k==max+1;
0096     disp('maximum number of iterations reached');
0097 end
0098 
0099 end
0100     
0101 
0102     
0103     
0104

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