function [xmin,fmin] = nr_brent(f, ax, bx, cx, tol)
%NR_BRENT - Brent's Method for Minimization.
%
% [ xmin, fmin ] = brent(f, ax, bx, cx, tol)
%
% Given a function f, and given a bracketing triplet of abscissas ax,
% bx and cx, and fx(b) is less than both f(ax) and f(cx), this routine
% isolates the minimum to a fractional precision of about TOL using
% Brent's method. The abscissa of the minimum is returned in xmin, and
% the minimum function value is returned in fmin.
%
% References:
%
% Brent, R.P. 1973, Algorithms for Minimization without Derivatives
% (Englewood Cliffs, NJ: Prentice Hall), Chapter 5
% Press, W.H. et al, 1992, Numerical recipes in C: the art of scientific
% computing (Cambridge University Press), Chapter 10.2
%
% Description of variables:
%
% a, b the minimum is bracketed between a and b
% x point with the least function value found so far
% w point with the second least function value found so far
% v previous value of w
% u point at which the function was evaluated most recently
% xm midpoint between a and b (the function is not evaluated there)
%
% Note: these points are not necessarily all distinct.
%
% e movement from best current value in the last iteration step
% etemp movement from the best current value in the second last
% iteration step
%
% General principles:
%
% - The parabola is fitted trough the points x, v and w
% - To be acceptable the parabolic step must
% (i) fall within the bounding interval (a,b)
% (ii) imply a movement from the best current value x that is
% *less* than half the movement of the *step before last*
% - The code never evaluates the function less than a distance tol1
% from a point already evaluated or from a known bracketing point
%
% Written by:
% --
% John L. Weatherwax 2004-12-11
%
% email: wax@alum.mit.edu
%
% Please send comments and especially bug reports to the
% above email address.
%
%-----
VERBOSE = 1; % Print steps taken for the solution.
CGOLD = (3-sqrt(5))/2; % golden ratio
ITMAX = 100; % max number of iterations
ZEPS = 1e-10; % absolute error tolerance
% initialization
e = 0;
if ax < cx
a = ax; b = cx;
else
a = cx; b = ax;
end;
v = bx; w = v; x = w;
fx = feval(f,x); fv = fx; fw = fv;
for iter = 1:ITMAX,
if( VERBOSE )
fprintf(1, 'k=%4d, |a-b|=%e\n', iter, abs(a-b));
end
xm = 0.5*(a+b);
tol1 = tol*abs(x) + ZEPS;
tol2 = 2*tol1;
% Stopping criterion: equivalent to: max(x-a, b-x) <= tol2
if abs(x-xm) <= tol2-0.5*(b-a)
xmin = x;
fmin = fx;
return
end
if abs(e) > tol1
%
% The second last move was sufficently large:
% let's construct the parabolic fit
%
r = (x-w)*(fx-fv);
q = (x-v)*(fx-fw);
p = (x-v)*q - (x-w)*r;
q = 2.0*(q-r);
if q > 0, p = -p; end
q = abs(q);
etemp = e;
e = d;
if abs(p) >= abs(0.5*q*etemp) | p <= q*(a-x) | p >= q*(b-x)
%
% The parabolic fit did not meet the reqirements above:
% use a golden section refinement instead.
%
% Print an explanation of what condition failed:
%
if( VERBOSE )
if abs(p) >= abs(0.5*q*etemp)
disp('abs(p) >= abs(0.5*q*etemp)')
end
if p <= q*(a-x)
disp('p <= q*(a-x)')
end
if p >= q*(b-x)
disp('p >= q*(b-x)')
end
end
%
% Set e in such a way, that a parabolic fit is possible in the
% next iteration.
%
if x >= xm
e = a-x;
else
e = b-x;
end
d = CGOLD*e;
if( VERBOSE )
disp('golden section step');
end
else
%
% the parabolic fit meets the above requirements: use iteration
%
d = p/q;
u = x+d;
if u-a < tol2 | b-u < tol2
%
% the parabolic fit is too close to the boundaries of the
% bracketing interval: put new point tol1 apart from the
% midpoint
%
d = tol1*sign(xm-x);
if( VERBOSE ) disp('too close to interval boundaries'); end
else
if( VERBOSE ) disp('parabolic step'); end
end
end
else
%
% The second last move was not sufficently large: use golden section
% set e in such a way, that a parabolic fit is possible in the next
% iteration.
%
if x >= xm
e = a-x;
else
e = b-x;
end
d = CGOLD*e;
if( VERBOSE )
disp('abs(e) > tol1');
disp('golden section step');
end
end
if abs(d) >= tol1
u = x+d;
else
%
% u is too close to x: put u farer apart
%
u = x + tol1*sign(d);
if( VERBOSE ) disp('u is too close to x'); end;
end
fu = feval(f,u); % finally evaluate the function
if fu <= fx
%
% u is better than x: x becomes new boundary
%
if u >= x
a = x;
else
b = x;
end
v = w; w = x; x = u;
fv = fw; fw = fx; fx = fu;
else
%
% x is better than u: u becomes new boundary
%
if u < x
a = u;
else
b = u;
end
if fu <= fw | w == x
v = w; w = u;
fv = fw; fw = fu;
elseif fu <= fv | v == x | v == w
v = u;
fv = fu;
end
end
end
error('too many iterations in brent');