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
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