An approximation x to the point where f attains a minimum on the interval (ax,bx) is determined.
the method used is a combination of golden section search and successive parabolic interpolation. convergence is never much slower than that for a fibonacci search. if f has a continuous second derivative which is positive at the minimum (which is not at ax or bx), then convergence is superlinear, and usually of the order of about 1.324.
the function f is never evaluated at two points closer together than epsabs(fmin) + (tol/3), where eps is approximately the square root of the relative machine precision. if f is a unimodal function and the computed values of f are always unimodal when separated by at least epsabs(x) + (tol/3), then fmin approximates the abcissa of the global minimum of f on the interval ax,bx with an error less than 3epsabs(fmin) + tol. if f is not unimodal, then fmin may approximate a local, but perhaps non-global, minimum to the same accuracy.
this function subprogram is a slightly modified version of the algol 60 procedure localmin given in richard brent, algorithms for minimization without derivatives, prentice - hall, inc. (1973).
[1] http://www.netlib.org/fmm/fmin.f
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(brent_class), | intent(inout) | :: | me | |||
real(kind=wp), | intent(in) | :: | ax |
left endpoint of initial interval |
||
real(kind=wp), | intent(in) | :: | bx |
right endpoint of initial interval |
||
real(kind=wp), | intent(in) | :: | tol |
desired length of the interval of uncertainty of the final result (>=0) |
abcissa approximating the point where f attains a minimum
function fmin(me,ax,bx,tol) result(xmin) implicit none class(brent_class),intent(inout) :: me real(wp),intent(in) :: ax !! left endpoint of initial interval real(wp),intent(in) :: bx !! right endpoint of initial interval real(wp),intent(in) :: tol !! desired length of the interval of !! uncertainty of the final result (>=0) real(wp) :: xmin !! abcissa approximating the point where !! f attains a minimum real(wp) :: a,b,d,e,xm,p,q,r,tol1,tol2,u,v,w real(wp) :: fu,fv,fw,fx,x real(wp) :: abs,sqrt,sign real(wp),parameter :: c = (three-sqrt(five))/two !! squared inverse of golden ratio real(wp),parameter :: half = 0.5_wp real(wp),parameter :: eps = sqrt(epsilon(one)) ! initialization a = ax b = bx v = a + c*(b - a) w = v x = v e = zero fx = me%f(x) fv = fx fw = fx do ! main loop starts here xm = half*(a + b) tol1 = eps*abs(x) + tol/three tol2 = two*tol1 ! check stopping criterion if (abs(x - xm) <= (tol2 - half*(b - a))) exit ! is golden-section necessary if (abs(e) <= tol1) then ! a golden-section step if (x >= xm) then e = a - x else e = b - x end if d = c*e else ! fit parabola r = (x - w)*(fx - fv) q = (x - v)*(fx - fw) p = (x - v)*q - (x - w)*r q = two*(q - r) if (q > zero) p = -p q = abs(q) r = e e = d ! is parabola acceptable if ((abs(p) >= abs(half*q*r)) .or. (p <= q*(a - x)) .or. (p >= q*(b - x))) then ! a golden-section step if (x >= xm) then e = a - x else e = b - x end if d = c*e else ! a parabolic interpolation step d = p/q u = x + d ! f must not be evaluated too close to ax or bx if (((u - a) < tol2) .or. ((b - u) < tol2)) d = sign(tol1, xm - x) end if end if ! f must not be evaluated too close to x if (abs(d) >= tol1) then u = x + d else u = x + sign(tol1, d) end if fu = me%f(u) ! update a, b, v, w, and x if (fu <= fx) then if (u >= x) a = x if (u < x) b = x v = w fv = fw w = x fw = fx x = u fx = fu else if (u < x) a = u if (u >= x) b = u if (fu <= fw .or. w == x) then v = w fv = fw w = u fw = fu else if (fu <= fv .or. v == x .or. v == w ) then v = u fv = fu end if end if end do ! end of main loop xmin = x end function fmin