real function fmin(ax,bx,f,tol) real ax,bx,f,tol c c an approximation x to the point where f attains a minimum on c the interval (ax,bx) is determined. c c c input.. c c ax left endpoint of initial interval c bx right endpoint of initial interval c f function subprogram which evaluates f(x) for any x c in the interval (ax,bx) c tol desired length of the interval of uncertainty of the final c result ( .ge. 0.0) c c c output.. c c fmin abcissa approximating the point where f attains a minimum c c c the method used is a combination of golden section search and c successive parabolic interpolation. convergence is never much slower c than that for a fibonacci search. if f has a continuous second c derivative which is positive at the minimum (which is not at ax or c bx), then convergence is superlinear, and usually of the order of c about 1.324.... c the function f is never evaluated at two points closer together c than eps*abs(fmin) + (tol/3), where eps is approximately the square c root of the relative machine precision. if f is a unimodal c function and the computed values of f are always unimodal when c separated by at least eps*abs(x) + (tol/3), then fmin approximates c the abcissa of the global minimum of f on the interval ax,bx with c an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, c then fmin may approximate a local, but perhaps non-global, minimum to c the same accuracy. c this function subprogram is a slightly modified version of the c algol 60 procedure localmin given in richard brent, algorithms for c minimization without derivatives, prentice - hall, inc. (1973). c c real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w real fu,fv,fw,fx,x c c c is the squared inverse of the golden ratio c c = 0.5*(3. - sqrt(5.0)) c c eps is approximately the square root of the relative machine c precision. c eps = 1.00 10 eps = eps/2.00 tol1 = 1.0 + eps if (tol1 .gt. 1.00) go to 10 eps = sqrt(eps) c c initialization c a = ax b = bx v = a + c*(b - a) w = v x = v e = 0.0 fx = f(x) fv = fx fw = fx c c main loop starts here c 20 xm = 0.5*(a + b) tol1 = eps*abs(x) + tol/3.0 tol2 = 2.0*tol1 c c check stopping criterion c if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 c c is golden-section necessary c if (abs(e) .le. tol1) go to 40 c c fit parabola c r = (x - w)*(fx - fv) q = (x - v)*(fx - fw) p = (x - v)*q - (x - w)*r q = 2.00*(q - r) if (q .gt. 0.0) p = -p q = abs(q) r = e e = d c c is parabola acceptable c 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 if (p .le. q*(a - x)) go to 40 if (p .ge. q*(b - x)) go to 40 c c a parabolic interpolation step c d = p/q u = x + d c c f must not be evaluated too close to ax or bx c if ((u - a) .lt. tol2) d = sign(tol1, xm - x) if ((b - u) .lt. tol2) d = sign(tol1, xm - x) go to 50 c c a golden-section step c 40 if (x .ge. xm) e = a - x if (x .lt. xm) e = b - x d = c*e c c f must not be evaluated too close to x c 50 if (abs(d) .ge. tol1) u = x + d if (abs(d) .lt. tol1) u = x + sign(tol1, d) fu = f(u) c c update a, b, v, w, and x c if (fu .gt. fx) go to 60 if (u .ge. x) a = x if (u .lt. x) b = x v = w fv = fw w = x fw = fx x = u fx = fu go to 20 60 if (u .lt. x) a = u if (u .ge. x) b = u if (fu .le. fw) go to 70 if (w .eq. x) go to 70 if (fu .le. fv) go to 80 if (v .eq. x) go to 80 if (v .eq. w) go to 80 go to 20 70 v = w fv = fw w = u fw = fu go to 20 80 v = u fv = fu go to 20 c c end of main loop c 90 fmin = x return end