Home > K25104 > NonlinearSystems > Brent.m

Brent

PURPOSE ^

implements Brent's algorithm to find a solution of f(x)=0;

SYNOPSIS ^

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

DESCRIPTION ^

 implements Brent's algorithm to find a solution of f(x)=0;
 Input arguments:
   f, function handle or vector of polynomial coefficients
   b, initial iterate
   c, initial contrapoint
   tol, tolerance
   max, maximum number of iterations
 Output arguments:
   x, solution
   k, number of iterations
   z, holds information which method was used for each iteration; 
     1 -> Binary search
     2 -> Inverse quadratic interpolation
     3 -> Linear interpolation

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ x,k,z ] = Brent( f,a,b,tol,max )
0002 % implements Brent's 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 %   c, initial contrapoint
0007 %   tol, tolerance
0008 %   max, maximum number of iterations
0009 % Output arguments:
0010 %   x, solution
0011 %   k, number of iterations
0012 %   z, holds information which method was used for each iteration;
0013 %     1 -> Binary search
0014 %     2 -> Inverse quadratic interpolation
0015 %     3 -> Linear interpolation
0016 
0017 % first check user inputs
0018 if tol<=0;
0019     error('tolerance must be positive');
0020 elseif max<=0;
0021     error('maximum number of iterations must be positive');
0022 elseif isa(f,'function_handle');
0023     % do nothing
0024 elseif isa(f,'double'); 
0025     [n,m]=size(f);  % find the size of f
0026     if n~=1 && m~=1; % then f is not a vector
0027         error(['input must be a function handle or vector holding ',...
0028             'polynomial coefficients']);
0029     else 
0030         v=poly2sym(f);  % this converts the vector of coefficients to
0031                         % an expression for the polynomial
0032         f=matlabFunction(v);    % this converts the expression to a
0033                                 % function handle
0034     end
0035 else
0036     error(['input must be a function handle or vector holding ',...
0037             'polynomial coefficients']);
0038 end
0039  
0040 delta=tol/2;        % set delta
0041 
0042 fa=feval(f,a);      % evaluate f(a)
0043 fb=feval(f,b);      % evaluate f(b)
0044 k=1;                % initialize iteration counter
0045 
0046 if fa==0;      % check in case the solution is at one of the endpoints
0047     x=a; return;
0048 elseif fb==0;
0049     x=b; return;
0050 elseif fa*fb>0;   % if f(a),f(b) have the same sign
0051         error('function values at endpoints have the same sign');
0052 end
0053 
0054 if abs(fb)<abs(fa);     % then we need to interchange a and b, so the
0055                         % iterate is the best approximation to x
0056     swap=a;             % hold a
0057     a=b; a=swap;        % swap a and b
0058     fa=feval(f,a);      % revaluate f(a)
0059     fb=feval(f,b);      % revaluate f(b)
0060 end
0061 
0062 c=a;            % we will need an extra variables to hold the previous
0063 d=a;            % iterates, in the first instance set these to a.
0064 fc=fa; fd=fa;       % store f(c), f(d)
0065 
0066 z=zeros(max,1);     % z will be a variable which holds information on
0067 z(1)=1;             % whether the previous iteration used binary
0068                     % search, linear interpolation or IQI.
0069                     
0070 
0071 
0072 % start of algorithm
0073 while abs(a-b)>tol && k<=max;
0074     
0075     if fa~=fb && fb~=fc && fc~=fa;  % then choose s according to
0076                                   % inverse quadratic interpolation
0077         s=a*fb*fc/((fa-fb)*(fa-fc))+fa*b*fc/((fb-fa)*(fb-fc))...
0078             +fa*fb*c/((fc-fa)*(fc-fb));
0079         temp=2;     % temporarily store this choice of interpolation
0080     else                          % if any of fa,fb,fc coincide,
0081                                   % choose s according to linear
0082                                   % interpolation
0083         s=b-fb*(b-a)/(fb-fa);
0084         temp=3;     % temporarily store this choice of interpolation
0085     end
0086     
0087     if or((s<(3*a+b)/4 & s<b),(s>(3*a+b)/4 & s>b)); % if s is not
0088                                                     % between (3*a+b)/4
0089                                                     % and b
0090         s=(a+b)/2; z(k+1)=1; % choose s according to binary search
0091                              % and set z=1
0092     elseif z(k)==1 && abs(b-c)<delta;
0093         s=(a+b)/2; z(k+1)=1; % choose s according to binary search
0094                              % and set z=1
0095     elseif z(k)~=1 && abs(c-d)<delta;
0096         s=(a+b)/2; z(k+1)=1; % choose s according to binary search
0097                              % and set z=1
0098     elseif z(k)==1 && abs(s-b)>=abs(b-a)/2;
0099         s=(a+b)/2; z(k+1)=1; % choose s according to binary search
0100                              % and set z=1
0101     elseif z(k)~=1 && abs(s-b)>=abs(b-d)/2;
0102         s=(a+b)/2; z(k+1)=1; % choose s according to binary search
0103                              % and set z=1
0104     else
0105         z(k+1)=temp;
0106     end
0107     
0108     fs=feval(f,s);     % evaluate f(s)
0109     if fs==0;          % check if we have found the solution
0110         x=s;
0111         return;
0112     end
0113      
0114     d=c;    % update d
0115     c=b;    % update c
0116     fa=feval(f,a);     % evaluate f(a)
0117     fb=feval(f,b);     % evaluate f(b)
0118     fd=feval(f,d);     % evaluate f(d)
0119     fc=feval(f,c);     % evaluate f(c)
0120     
0121     if fa*fs<0;   % if f(a) and f(s) have different signs
0122         b=s;        % update b
0123         fb=feval(f,b);     % update f(b)
0124     else 
0125         a=s;        % update a
0126         fa=feval(f,a);     % update f(a)
0127     end
0128     
0129     if abs(fa)<abs(fb);     % then interchange a and b, so the iterate
0130                             % is the best approximation to x
0131         swap=a;             % hold a
0132         a=b; b=swap;        % swap b and a
0133         swap=fa;            % hold f(a)
0134         fa=fb;  fb=swap;    % swap f(a) and f(b)
0135     end
0136  
0137     if fb==0;      % check if we have found the solution
0138         x=a;
0139     end
0140     
0141     k=k+1;  % increment k
0142 end
0143 
0144 x=b;        % set x to the most recent iterate
0145 z=z(1:k);   % shorten z
0146 
0147 if k==max+1;
0148     disp('maximum number of iterations reached');
0149 end
0150     
0151 end

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