function [xmin,fmin] = nr_dbrent(f, fp, ax, bx, cx, tol) %NR_DBRENT - Brent's Method for Minimization. % % [ xmin, fmin ] = dbrent(f, fp, ax, bx, cx, tol) % % Given a function f and its derivative function df, and given a % bracketing triplet of abscissas ax, bx and cx [ such that bx is % between ax 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 that uses derivatives. The abscissa of the % minimum is returned in xmin, and the minimum function value is returned % in fmin. % % References: % % Press, W.H. et al, 1992, Numerical recipes in FORTRAN: 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 % VERBOSE = 1; % Print steps taken for the solution. 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; dx = feval(fp,x); dv = dx; dw = dx; 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 ) d1 = 2.0*(b-a); d2 = d1; if( dw ~= dx ) d1 = (w-x)*dx/(dx-dw); % Secant method with one point. end if( dv ~= dx ) d2 = (v-x)*dx/(dx-dv); % And the other. end u1 = x+d1; u2 = x+d2; ok1 = (a-u1)*(u1-b) > 0.0 & dx*d1 <= 0.0; ok2 = (a-u2)*(u2-b) > 0.0 & dx*d2 <= 0.0; olde=e; e=d; if( ok1 | ok2 ) if( ok1 & ok2 ) if( abs(d1) < abs(d2) ) d = d1; else d = d2; end else if( ok1 ) d = d1; else d = d2; end end if( abs(d) <= abs(0.5*olde) ) u=x+d; if( u-a < tol2 | b-u < tol2 ) d = SIGN(tol1,xm-x); end else if( dx >= 0.0 ) e = a-x; else e = b-x; end d = 0.5*e; end else if( dx >= 0.0 ) e = a-x; else e = b-x; end d = 0.5*e; end else if( dx >= 0.0 ) e = a-x; else e = b-x; end d = 0.5*e; end % End if( abs(e) > tol1 ) if( abs(d) >= tol1 ) u = x+d; fu = feval(f,u); else % % u is too close to x: put u farther apart % u = x + SIGN(tol1,d); fu = feval(f,u); if( fu > fx ) xmin = x; fx = fu; end end du = feval(fp,u); if( fu <= fx ) % % u is better than x: x becomes new boundary % if( u >= x ) a = x; else b = x; end v = w; fv = fw; dv = dw; w = x; fw = fx; dw = dx; x = u; fx = fu; dx = du; 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; fv = fw; dv = dw; w = u; fw = fu; dw = du; else if( fu <= fv | v == x | v == w ) v = u; fv = fu; dv = du; end end end end % End for loop. error('Too many iterations in dbrent.');