home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / FZERO.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.2 KB  |  107 lines
  1. function b = fzero(FunFcn,x,tol,trace)
  2. %FZERO    Find a zero of a function of one variable. 
  3. %    FZERO(F,X) finds a zero of f(x).  F is a string containing the
  4. %    name of a real-valued function of a single real variable.   X is
  5. %    a starting guess.  The value returned is near a point where F
  6. %    changes sign.  For example, FZERO('sin',3) is pi.  Note the
  7. %    quotes around sin.  Ordinarily, functions are defined in M-files.
  8. %
  9. %    An optional third argument sets the relative tolerance for the
  10. %    convergence test.   The presence of an nonzero optional fourth
  11. %    argument triggers a printing trace of the steps.
  12.  
  13. %    C.B. Moler 1-19-86
  14. %    Revised CBM 3-25-87, LS 12-01-88.
  15. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  16.  
  17. %  This algorithm was originated by T. Dekker.  An Algol 60 version,
  18. %  with some improvements, is given by Richard Brent in "Algorithms for
  19. %  Minimization Without Derivatives", Prentice-Hall, 1973.  A Fortran
  20. %  version is in Forsythe, Malcolm and Moler, "Computer Methods
  21. %  for Mathematical Computations", Prentice-Hall, 1976.
  22.  
  23. % Initialization
  24.  
  25. if nargin < 3, trace = 0; tol = eps; end
  26. if nargin == 3, trace = 0; end
  27. if trace, clc, end
  28. if (length(x) > 1) | (~finite(x))
  29.    error('Second argument must be a finite scalar.')
  30. end
  31. if x ~= 0, dx = x/20;
  32. else, dx = 1/20;
  33. end
  34. a = x - dx;  fa = feval(FunFcn,a);
  35. if trace, home, init = [a fa], end
  36. b = x + dx;  fb = feval(FunFcn,b);
  37. if trace, home, init = [b fb], end
  38.  
  39. % Find change of sign.
  40.  
  41. while (fa > 0) == (fb > 0)
  42.    dx = 2*dx;
  43.    a = x - dx;  fa = feval(FunFcn,a);
  44.    if trace, home, sign = [a fa], end
  45.    if (fa > 0) ~= (fb > 0), break, end
  46.    b = x + dx;  fb = feval(FunFcn,b);
  47.    if trace, home, sign = [b fb], end
  48. end
  49.  
  50. fc = fb;
  51. % Main loop, exit from middle of the loop
  52. while fb ~= 0
  53.    % Insure that b is the best result so far, a is the previous
  54.    % value of b, and c is on the opposite of the zero from b.
  55.    if (fb > 0) == (fc > 0)
  56.       c = a;  fc = fa;
  57.       d = b - a;  e = d;
  58.    end
  59.    if abs(fc) < abs(fb)
  60.       a = b;    b = c;    c = a;
  61.       fa = fb;  fb = fc;  fc = fa;
  62.    end
  63.  
  64.    % Convergence test and possible exit
  65.    m = 0.5*(c - b);
  66.    toler = 2.0*tol*max(abs(b),1.0);
  67.    if (abs(m) <= toler) + (fb == 0.0), break, end
  68.  
  69.    % Choose bisection or interpolation
  70.    if (abs(e) < toler) + (abs(fa) <= abs(fb))
  71.    % Bisection
  72.       d = m;  e = m;
  73.    else
  74.    % Interpolation
  75.       s = fb/fa;
  76.       if (a == c)
  77.       % Linear interpolation
  78.          p = 2.0*m*s;
  79.          q = 1.0 - s;
  80.       else
  81.       % Inverse quadratic interpolation
  82.          q = fa/fc;
  83.          r = fb/fc;
  84.          p = s*(2.0*m*q*(q - r) - (b - a)*(r - 1.0));
  85.          q = (q - 1.0)*(r - 1.0)*(s - 1.0);
  86.       end;
  87.       if p > 0, q = -q; else p = -p; end;
  88.       % Is interpolated point acceptable
  89.       if (2.0*p < 3.0*m*q - abs(toler*q)) * (p < abs(0.5*e*q))
  90.          e = d;  d = p/q;
  91.       else
  92.          d = m;  e = m;
  93.       end;
  94.    end % Interpolation
  95.  
  96.    % Next point
  97.    a = b;
  98.    fa = fb;
  99.    if abs(d) > toler, b = b + d;
  100.    else if b > c, b = b - toler;
  101.         else b = b + toler;
  102.         end
  103.    end
  104.    fb = feval(FunFcn,b);
  105.    if trace, home, step = [b fb], end
  106. end % Main loop
  107.