function [ ax, bx, cx, fa, fb, fc ] = nr_mnbrak(ax,bx,func) % NR_MNBRAK - returns points which bracket a minimum of the function % % Given a function func, and given distinct initial points ax and bx, % this routine searches in the downhill direction (defined by the function % as evaluated at the initial points) and returns new points ax, bx, cx, % which bracket a minimum of the function. Also returned are the % function values at the three points. % Note: The use of the same name ax, and bx on both the input and the % output of this function requires that both values are COPIED into % the functions workspace. In general elements listed as inputs only are % simply read and not copied. % % 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. % %----- GOLD = 1.618034; G_LIMIT = 100.0; TINY = 1.d-20; fa = feval(func,ax); fb = feval(func,bx); if( fb > fa ) dum = ax; ax = bx; bx = dum; dum = fb; fb = fa; fa = dum; end; cx = bx+GOLD*(bx-ax); fc = feval(func,cx); while( fb >= fc ) r=(bx-ax)*(fb-fc); q=(bx-cx)*(fb-fa); % u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*(sign(q-r)+eps)*max(abs(q-r),TINY)); u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*SIGN(max(abs(q-r),TINY),q-r)); ulim = bx+G_LIMIT*(cx-bx); if( (bx-u)*(u-cx) > 0.0 ) fu = feval(func,u); if( fu < fc ) ax = bx; fa = fb; bx = u; fb = fu; return; elseif( fu > fb ) cx = u; fc = fu; return; end; u=cx+GOLD*(cx-bx); fu = feval(func,u); elseif( (cx-u)*(u-ulim) > 0.0 ) fu = feval(func,u); if( fu < fc ) bx = cx; cx=u; u=cx+GOLD*(cx-bx); fb = fc; fc = fu; fu = feval(func,u); end; elseif( (u-ulim)*(ulim-cx) >= 0.0 ) u=ulim; fu = feval(func,u); else u=cx+GOLD*(cx-bx); fu = feval(func,u); end; ax = bx; bx = cx; cx = u; fa = fb; fb = fc; fc = fu; end;