/******************************************************* * Find minimum of a real function Y=F(X) * * ---------------------------------------------------- * * REFERENCE: "Numerical Recipes, The Art of Scientific * * Computing By W.H. Press, B.P. Flannery, * * S.A. Teukolsky and W.T. Vetterling, * * Cambridge University Press, 1986" * * [BIBLI 08]. * * ---------------------------------------------------- * * Sample run: * * * * Find a minimum of function X*SIN(X)-2*COS(X) * * between X=4 and X=6. * * * * X1=4 X2=5 X3=6 TOL=1e-6 * * * * Function minimum is: -5.534529 * * for X= 5.232937 * * * *Press any key to continue * * * * ---------------------------------------------------- * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * *******************************************************/ #include #include #define ITMAX 100 #define CGOLD 0.3819660 #define ZEPS 1.0E-10 //Maximum allowed number of iterations; golden ration; and a small //number which protects against trying to achieve fractional accuracy //for a minimum that happens to be exactly zero. double x1,x2,x3,xmini,ymini,tol; //Function to be analyzed double F(double x) { return x*sin(x)-2.0*cos(x); } double MIN(double a, double b) { if (a<=b) return a; else return b; } double MAX(double a, double b) { if (a>=b) return a; else return b; } double Sign(double a, double b) { if (b>=0) return fabs(a); else return -fabs(a); } double BRENT(double ax, double bx, double cx, double tol, double *xmin) { // Labels: e1,e2,e3 double a,b,d,e,fu,fv,fx,fw,u,v,w,x; double etemp,p,q,r,tol1,tol2,xm; int iter; /*Given a function F, and 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 in XMIN, and the minimum function value is returned as BRENT, the returned function value.*/ a=MIN(ax,cx); b=MAX(ax,cx); v=bx; w=v; x=v; e=0; fx=F(x); fv=fx; fw=fx; d=0; for (iter=1; iter<=ITMAX; iter++) { // main loop xm=0.5*(a+b); tol1=tol*fabs(x)+ZEPS; tol2=2.0*tol1; if (fabs(x-xm)<=tol2-0.5*(b-a)) goto e3; //Test for done here 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); //bug corrected 07/24/2006 (0.2 instead of 2.0) if (q>0) p=-p; q=fabs(q); etemp=e; e=d; if (fabs(p)>=fabs(0.5*q*etemp) || p<=q*(a-x) || p>=q*(b-x)) goto e1; // The above conditions determine the acceptability of the // parabolic fit. Here it is o.k.: d=p/q; u=x+d; if (u-a=xm) e=a-x; else e=b-x; d=CGOLD*e; e2: if (fabs(d)>=tol1) u=x+d; else u=x+Sign(tol1,d); fu=F(u); // This the one function evaluation per iteration if (fu<=fx) { if (u>=x) a=x; else b=x; v=w; fv=fw; w=x; fw=fx; x=u; fx=fu; } else { if (u