/******************************************************* * Bracketing a minimum of a real function Y=F(X) * * using MNBRAK subroutine * * ---------------------------------------------------- * * 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 3 points of function X*SIN(X)-2*COS(X) * * bracketing a minimum, given two initial points. * * * * X1=4.0 X2=5.0 * * * * The three points are: * * * * X1= 4.000000 X2= 5.000000 X3= 6.618034 * * * * Corresponding function values: * * * * F1= -1.719923 F2= -5.361946 F3= 0.285940 * * * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * *******************************************************/ #include #include /*The first parameter is the default ratio by which successive intervals are magnified; the second is the maximum magnification allowed for a parabolic-fit step. */ #define GOLD 1.618034 #define GLIMIT 100 #define TINY 1e-20 double x1, x2, x3, f1, f2, f3; // Function to be analyzed double FUNC(double x) { return x*sin(x)-2*cos(x); } 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); } void MNBRAK(double ax, double bx, double *cx, double *fa, double *fb, double *fc) { /*Given a function FUNC(X), 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, fa, fb and fc. */ // LABEL: e1 double dum,fu,q,r,u,ulim; *fa=FUNC(ax); *fb=FUNC(bx); if (*fb > *fa) { dum=ax; ax=bx; bx=dum; dum=*fb; *fb=*fa; *fa=dum; } *cx=bx+GOLD*(bx-ax); *fc=FUNC(*cx); e1:if (*fb>=*fc) { r=(bx-ax)*(*fb-*fc); q=(bx-*cx)*(*fb-*fa); u=bx-((bx-*cx)*q-(bx-ax)*r)/(2.*Sign(Max(fabs(q-r),TINY),q-r)); ulim=bx+GLIMIT*(*cx-bx); if ((bx-u)*(u-*cx)>0) { fu=FUNC(u); if (fu < *fc) { ax=bx; *fa=*fb; bx=u; *fb=fu; goto e1; } else if (fu > *fb) { *cx=u; *fc=fu; goto e1; } u=*cx+GOLD*(*cx-bx); fu=FUNC(u); } else if ((*cx-u)*(u-ulim) > 0) { fu=FUNC(u); if (fu < *fc) { bx=*cx; *cx=u; u=*cx+GOLD*(*cx-bx); *fb=*fc; *fc=fu; fu=FUNC(u); } } else if ((u-ulim)*(ulim-*cx)>=0) { u=ulim; fu=FUNC(u); } else { u=*cx+GOLD*(*cx-bx); fu=FUNC(u); } ax=bx; bx=*cx; *cx=u; *fa=*fb; *fb=*fc; *fc=fu; goto e1; } } void main() { x1=4.0; x2=5.0; MNBRAK(x1,x2,&x3,&f1,&f2,&f3); printf("\n The three points are:\n\n"); printf(" X1=%10.6f X2=%10.6f X3=%10.6f\n",x1,x2,x3); printf("\n Corresponding function values:\n\n"); printf(" F1=%10.6f F2=%10.6f F3=%10.6f\n\n\n",f1,f2,f3); } // end of file mnbrak.cpp