home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / FMIN.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.2 KB  |  141 lines
  1. function [xf,options] = fmin(funfcn,ax,bx,options,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
  2. %FMIN    Minimize a function of one variable.
  3. %    FMIN('F',x1,x2) attempts to return a value of x which is a local 
  4. %    minimizer of F(x) in the interval x1 < x < x2.  'F' is a string 
  5. %    containing the name of the objective function to be minimized.
  6. %
  7. %    FMIN('F',x1,x2,OPTIONS) uses a vector of control parameters.
  8. %    If OPTIONS(1) is nonzero, intermediate steps in the solution are
  9. %    displayed; the default is OPTIONS(1) = 0.  OPTIONS(2) is the termination
  10. %    tolerance for x; the default is 1.e-4.  OPTIONS(14) is the maximum
  11. %    number of steps; the default is OPTIONS(14) = 500.  The other components
  12. %    of OPTIONS are not used as input control parameters by FMIN.  For more
  13. %    information, see FOPTIONS.
  14. %
  15. %    FMIN('F',x1,x2,OPTIONS,P1,P2,...) provides for up to 10 additional
  16. %    arguments which are passed to the objective function, F(X,P1,P2,...)
  17. %
  18. %    Example:
  19. %        fmin('cos',3,4) computes pi to a few decimal places.
  20. %        fmin('cos',3,4,[1,1.e-12]) displays the steps taken
  21. %        to compute pi to about 12 decimal places.  
  22.  
  23. %    Reference: "Computer Methods for Mathematical Computations",
  24. %    Forsythe, Malcolm, and Moler, Prentice-Hall, 1976.
  25.  
  26. %    Revised 10-14-88 JNL, 7-9-90 ACWG, 1-17-92, 6-15-92 CBM.
  27. %    Original coding by Duane Hanselman, University of Maine.
  28. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  29.  
  30. % initialization
  31. if nargin<4, options = []; end
  32. options = foptions(options);
  33. print = options(1);
  34. tol = options(2);
  35.  
  36. evalstr = [funfcn];
  37. if ~any(funfcn<48)
  38.     evalstr=[evalstr, '(x'];
  39.     for i=1:nargin - 4
  40.         evalstr = [evalstr,',P',int2str(i)];
  41.     end
  42.     evalstr = [evalstr, ')'];
  43. end
  44.  
  45. if (~options(14))
  46.     options(14) = 500; 
  47. end
  48.  
  49. num = 1;
  50. seps = sqrt(eps);
  51. c = 0.5*(3.0 - sqrt(5.0));
  52. a = ax; b = bx;
  53. v = a + c*(b-a);
  54. w = v; xf = v;
  55. d = 0.0; e = 0.0;
  56. x= xf; fx = eval(evalstr);  
  57. if print, clc, fmin_data = [1 xf fx], end
  58. fv = fx; fw = fx;
  59. xm = 0.5*(a+b);
  60. tol1 = seps*abs(xf) + tol/3.0;   tol2 = 2.0*tol1;
  61.  
  62. % Main loop
  63.  
  64. options(10) = 0; 
  65. while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) )
  66.  
  67.     num = num+1;
  68.     gs = 1;
  69.  
  70.     % Is a parabolic fit possible
  71.     if abs(e) > tol1
  72.         % Yes, so fit parabola
  73.         gs = 0;
  74.         r = (xf-w)*(fx-fv);
  75.         q = (xf-v)*(fx-fw);
  76.         p = (xf-v)*q-(xf-w)*r;
  77.         q = 2.0*(q-r);
  78.         if q > 0.0,  p = -p; end
  79.         q = abs(q);
  80.         r = e;  e = d;
  81.  
  82.         % Is the parabola acceptable
  83.         if ( (abs(p)<abs(0.5*q*r)) & (p>q*(a-xf)) & (p<q*(b-xf)) )
  84.  
  85.             % Yes, parabolic interpolation step
  86.             d = p/q;
  87.             x = xf+d;
  88.             step = '   num        xf        fx       parabolic';
  89.      
  90.             % f must not be evaluated too close to ax or bx
  91.             if ((x-a) < tol2) | ((b-x) < tol2)
  92.                 si = sign(xm-xf) + ((xm-xf) == 0);
  93.                 d = tol1*si;
  94.             end
  95.         else
  96.             % Not acceptable, must do a golden section step
  97.             gs=1;
  98.         end
  99.     end
  100.     if gs
  101.         % A golden-section step is required
  102.         if xf >= xm, e = a-xf;    else, e = b-xf;  end
  103.             d = c*e;
  104.             step = '   num        xf        fx       golden   ';
  105.     end
  106.  
  107.     % The function must not be evaluated too close to xf
  108.     si = sign(d) + (d == 0);
  109.     x = xf + si * max( abs(d), tol1 );
  110.     fu = eval(evalstr);  
  111.     if print, clc, fmin_data = [num x fu], disp(step), end
  112.  
  113.     % Update a, b, v, w, x, xm, tol1, tol2
  114.     if fu <= fx
  115.         if x >= xf, a = xf; else, b = xf; end
  116.         v = w; fv = fw;
  117.         w = xf; fw = fx;
  118.         xf = x; fx = fu;
  119.     else % fu > fx
  120.         if x < xf, a = x; else,b = x; end
  121.         if ( (fu <= fw) | (w == xf) )
  122.             v = w; fv = fw;
  123.             w = x; fw = fu;
  124.         elseif ( (fu <= fv) | (v == xf) | (v == w) )
  125.             v = x; fv = fu;
  126.         end
  127.     end
  128.     xm = 0.5*(a+b);
  129.     tol1 = seps*abs(xf) + tol/3.0; tol2 = 2.0*tol1;
  130.     if num > options(14)
  131.         if options(1)>=0
  132.         disp('Warning: Maximum number of iterations has been exceeded');
  133.         disp('       - increase options(14) for more iterations.')
  134.         options(10)=num;
  135.         return
  136.         end
  137.     end
  138. end
  139. options(8) = eval(evalstr);
  140. options(10)=num;
  141.