404 Chapter 10.Minimization or Maximization of Functions useless,the method will approximately alternate between parabolic steps and golden sections,converging in due course by virtue of the latter.The reason for comparing to the step before last seems essentially heuristic:Experience shows that it is better not to"punish"the algorithm for a single bad step if it can make it up on the next one. Another principle exemplified in the code is never to evaluate the function less than a distance tol from a point already evaluated (or from a known bracketing point).The reason is that,as we saw in equation (10.1.2),there is simply no information content in doing so:the function will differ from the value already evaluated only by an amount of order the roundofferror.Therefore in the code below you will find several tests and modifications of a potential new point,imposing this 8 restriction.This restriction also interacts subtly with the test for "doneness,"which the method takes into account. A typical ending configuration for Brent's method is that a and b are 2 x xx tol apart,with(the best abscissa)at the midpoint of a and b,and therefore fractionally accurate to±tol Indulge us a final reminder that tol should generally be no smaller than the 3IE square root of your machine's floating-point precision. 6 #include <math.h> #include "nrutil.h" #define ITMAX 100 Americ computer, #define CGOLD 0.3819660 ART #define ZEPS 1.0e-10 Here ITMAX is the maximum allowed number of iterations:CGOLD is the golden ratio;ZEPS is 9 Programs a small number that protects against trying to achieve fractional accuracy for a minimum that happens to be exactly zero. #define SHFT(a,b,c,d)(a)=(b);(b)=(c);(c)=(d); float brent(float ax,float bx,float cx,float (*f)(float),float tol, float *xmin) to dir Given a function f,and given a bracketing triplet of abscissas ax,bx,cx (such that bx is between ax and cx.and f(bx)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 SCIENTIFIC COMPUTING(ISBN the minimum is returned as xmin,and the minimum function value is returned as brent,the 19891892 returned function value. int iter; 10-621 float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; float e=0.0; This will be the distance moved on idge.org Fuurggoglrion Numerical Recipes -43108 the step before last. a=(ax cx ax cx); a and b must be in ascending order b=(ax>cx？ax:cx); but input abscissas need not be. x=w-v=bx; Initializations... (outside fw=fv=fx=(*f)(x); Software. for (iter=1;iter<=ITMAX;iter++) Main program loop. Xm=0.5*(a+b); tol2=2.0*(tol1=tol*fabs(x)+ZEPS); if (fabs(x-xm)<=(tol2-0.5*(b-a))){ Test for done here. *xmin=x; return fx; if (fabs(e)>tol1){ Construct a trial parabolic fit. r=(x-w)*(fx-fv); q=(x-v)*(fx-fw); p=(x-V)*q-(x-W)*r; q=2.0*(q-r); 1f(q>0.0)p=-P q=fabs(q);404 Chapter 10. Minimization or Maximization of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). useless, the method will approximately alternate between parabolic steps and golden sections, converging in due course by virtue of the latter. The reason for comparing to the step before last seems essentially heuristic: Experience shows that it is better not to “punish” the algorithm for a single bad step if it can make it up on the next one. Another principle exemplified in the code is never to evaluate the function less than a distance tol from a point already evaluated (or from a known bracketing point). The reason is that, as we saw in equation (10.1.2), there is simply no information content in doing so: the function will differ from the value already evaluated only by an amount of order the roundoff error. Therefore in the code below you will find several tests and modifications of a potential new point, imposing this restriction. This restriction also interacts subtly with the test for “doneness,” which the method takes into account. A typical ending configuration for Brent’s method is that a and b are 2×x×tol apart, with x (the best abscissa) at the midpoint of a and b, and therefore fractionally accurate to ±tol. Indulge us a final reminder that tol should generally be no smaller than the square root of your machine’s floating-point precision. #include <math.h> #include "nrutil.h" #define ITMAX 100 #define CGOLD 0.3819660 #define ZEPS 1.0e-10 Here ITMAX is the maximum allowed number of iterations; CGOLD is the golden ratio; ZEPS is a small number that protects against trying to achieve fractional accuracy for a minimum that happens to be exactly zero. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin) Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is between ax and cx, and f(bx) 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 as xmin, and the minimum function value is returned as brent, the returned function value. { int iter; float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; float e=0.0; This will be the distance moved on the step before last. a=(ax < cx ? ax : cx); a and b must be in ascending order, b=(ax > cx ? ax : cx); but input abscissas need not be. x=w=v=bx; Initializations... fw=fv=fx=(*f)(x); for (iter=1;iter<=ITMAX;iter++) { Main program loop. xm=0.5*(a+b); tol2=2.0*(tol1=tol*fabs(x)+ZEPS); if (fabs(x-xm) <= (tol2-0.5*(b-a))) { Test for done here. *xmin=x; return fx; } if (fabs(e) > tol1) { Construct a trial 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.0) p = -p; q=fabs(q);