fzero_old.m (9679B)
1 function b = fzero_old(FunFcn,x,tol,trace,varargin) 2 %FZERO Find zero of function of one variable. 3 % The Matlab 5.2 Version (replaced in version 5.3) 4 % FZERO(F,X) tries to find a zero of F. F is a string containing 5 % the name of a real-valued function of a single real variable. 6 % The value returned is near a point where F changes sign, or NaN 7 % if the search fails. 8 % 9 % FZERO(F,X), where X is a vector of length 2, assumes X is an 10 % interval where the sign of F(X(1)) differs from the sign of F(X(2)). 11 % An error occurs if this is not true. Calling FZERO with an interval 12 % guarantees FZERO will return a value near a point where F changes 13 % sign. 14 % 15 % FZERO(F,X), where X is a scalar value, uses X as a starting guess. 16 % FZERO looks for an interval containing a sign change for F and 17 % containing X. If no such interval is found, NaN is returned. 18 % In this case, the search terminates when the search interval 19 % is expanded until an Inf, NaN, or complex value is found. 20 % 21 % FZERO(F,X,TOL) sets the relative tolerance for the convergence test. 22 % FZERO(F,X,TOL,TRACE) displays information at each iteration when 23 % TRACE is nonzero. 24 % FZERO(F,X,TOL,TRACE,P1,P2,...) allows for additional arguments 25 % which are passed to the function, F(X,P1,P2,...). Pass an empty 26 % matrix for TOL or TRACE to use the default value. 27 % 28 % Examples 29 % fzero('sin', 3) returns pi. Note the quotes around sin. 30 % Ordinarily, functions are defined in M-files. 31 % fzero('abs(x)+1', 1) returns NaN since this function does 32 % not change sign anywhere on the real axis (and does not have 33 % a zero as well). 34 % fzero('sin', 3, [], 1) returns pi, uses the default tolerance, 35 % and displays iteration information. 36 % 37 % See also ROOTS. 38 39 % Copyright (c) 1984-98 by The MathWorks, Inc. 40 % $Revision: 5.12 $ $Date: 1997/11/21 23:30:46 $ 41 42 % This algorithm was originated by T. Dekker. An Algol 60 version, 43 % with some improvements, is given by Richard Brent in "Algorithms for 44 % Minimization Without Derivatives", Prentice-Hall, 1973. A Fortran 45 % version is in Forsythe, Malcolm and Moler, "Computer Methods 46 % for Mathematical Computations", Prentice-Hall, 1976. 47 48 % Initialization 49 50 if nargin < 3 | isempty(tol), tol = eps; end 51 if nargin < 4 | isempty(trace), trace = 0; end 52 if trace 53 header = ' Func evals x f(x) Procedure'; 54 step=' '; 55 count = 0; 56 end 57 if (~isfinite(x)) 58 error('Second argument must be finite.') 59 end 60 61 % Convert to inline function as needed. 62 FunFcn = fcnchk(FunFcn,length(varargin)); 63 64 % Interval input 65 if (length(x) == 2) 66 a = x(1); 67 b = x(2); 68 fa = feval(FunFcn,a,varargin{:}); 69 fb = feval(FunFcn,b,varargin{:}); 70 if any(~isfinite([fa fb])) | any(~isreal([fa fb])) 71 error('Function values at interval endpoints must be finite and real.') 72 end 73 if trace 74 disp(header) 75 data = [a fa]; step=' initial'; 76 count = count + 1; 77 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 78 data = [b fb]; step = ' initial'; 79 count = count + 1; 80 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 81 end 82 if (fa > 0) == (fb > 0) 83 error('The function values at the interval endpoints must differ in sign.') 84 end 85 86 % Starting guess scalar input 87 elseif (length(x) == 1) 88 fx = feval(FunFcn,x,varargin{:}); 89 if fx == 0 90 b = x; 91 return 92 elseif ~isfinite(fx) | ~isreal(fx) 93 error('Function value at starting guess must be finite and real.'); 94 end 95 if trace 96 disp(header) 97 data = [x fx]; step=' initial'; 98 count = count + 1; 99 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 100 end 101 if x ~= 0, 102 dx = x/50; 103 else, 104 dx = 1/50; 105 end 106 107 % Find change of sign. 108 twosqrt = sqrt(2); 109 a = x; fa = fx; b = x; fb = fx; 110 111 while (fa > 0) == (fb > 0) 112 dx = twosqrt*dx; 113 a = x - dx; fa = feval(FunFcn,a,varargin{:}); 114 if ~isfinite(fa) | ~isreal(fa) 115 disperr(a,fa); 116 b = NaN; 117 return 118 end 119 if trace 120 data = [a fa]; step=' search'; 121 count = count + 1; 122 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 123 end 124 if (fa > 0) ~= (fb > 0) 125 break 126 end 127 b = x + dx; fb = feval(FunFcn,b,varargin{:}); 128 if ~isfinite(fb) | ~isreal(fb) 129 disperr(b,fb); 130 b = NaN; 131 return 132 end 133 if trace 134 data = [b fb]; step=' search'; 135 count = count + 1; 136 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 137 end 138 end % while 139 if trace 140 disp(' ') 141 disp([' Looking for a zero in the interval [', ... 142 num2str(a) , ', ', num2str(b), ']']); 143 disp(' ') 144 end 145 else 146 error('Second argument must be of length 1 or 2.'); 147 end % if (length(x) == 2 148 149 fc = fb; 150 % Main loop, exit from middle of the loop 151 while fb ~= 0 152 % Insure that b is the best result so far, a is the previous 153 % value of b, and c is on the opposite of the zero from b. 154 if (fb > 0) == (fc > 0) 155 c = a; fc = fa; 156 d = b - a; e = d; 157 end 158 if abs(fc) < abs(fb) 159 a = b; b = c; c = a; 160 fa = fb; fb = fc; fc = fa; 161 end 162 163 % Convergence test and possible exit 164 m = 0.5*(c - b); 165 toler = 2.0*tol*max(abs(b),1.0); 166 if (abs(m) <= toler) + (fb == 0.0), break, end 167 168 % Choose bisection or interpolation 169 if (abs(e) < toler) + (abs(fa) <= abs(fb)) 170 % Bisection 171 d = m; e = m; 172 step=' bisection'; 173 else 174 % Interpolation 175 s = fb/fa; 176 if (a == c) 177 % Linear interpolation 178 p = 2.0*m*s; 179 q = 1.0 - s; 180 else 181 % Inverse quadratic interpolation 182 q = fa/fc; 183 r = fb/fc; 184 p = s*(2.0*m*q*(q - r) - (b - a)*(r - 1.0)); 185 q = (q - 1.0)*(r - 1.0)*(s - 1.0); 186 end; 187 if p > 0, q = -q; else p = -p; end; 188 % Is interpolated point acceptable 189 if (2.0*p < 3.0*m*q - abs(toler*q)) * (p < abs(0.5*e*q)) 190 e = d; d = p/q; 191 step=' interpolation'; 192 else 193 d = m; e = m; 194 step=' bisection'; 195 end; 196 end % Interpolation 197 198 % Next point 199 a = b; 200 fa = fb; 201 if abs(d) > toler, b = b + d; 202 else if b > c, b = b - toler; 203 else b = b + toler; 204 end 205 end 206 fb = feval(FunFcn,b,varargin{:}); 207 if trace 208 data = [b fb]; 209 count = count + 1; 210 disp([sprintf('%5.0f %13.6g %13.6g ',count, data), step]) 211 end 212 end % Main loop 213 214 %------------------------------------------------------------------ 215 216 function disperr(y, fy) 217 %DISPERR Display an appropriate error message when FY is Inf, 218 % NaN, or complex. Assumes Y is the value and FY is the function 219 % value at Y. If FY is neither Inf, NaN, or complex, it generates 220 % an error message. 221 222 if ~isfinite(fy) % NaN or Inf detected 223 disp('NaN or Inf function value encountered during '); 224 disp(' search for an interval containing a sign change.'); 225 disp(['Function value at ', num2str(y),' is ',num2str(fy)]); 226 disp('Aborting since no such interval was found.') 227 disp('Check function or try again with a different starting value.') 228 elseif ~isreal(fy) % Complex value detected 229 disp('Complex function value encountered during '); 230 disp(' search for an interval containing a sign change.'); 231 disp(['Function value at ', num2str(y),' is ',num2str(fy)]); 232 disp('Aborting since no such interval was found.') 233 disp('Check function or try again with a different starting value.') 234 else 235 error('Disperr called with invalid argument.') 236 end