'******************************************************** '* 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" [8]. * '* ---------------------------------------------------- * '* 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 * '* * '* ---------------------------------------------------- * '* BASIC version by J-P Moreau. * '* (www.jpmoreau.fr) * '******************************************************** defdbl a-h,o-z ITMAX=100 CGOLD=0.3819660 ZEPS=1E-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. ax=4.0 : bx=5.0 : cx=6.0 tol=1e-6 gosub 2000 'call Brent method cls print print using " Function minimum is ##.######"; fx print print using " for X= ##.######"; x print : print end 'Function to be analyzed 1000 F=xx*sin(xx)-2*cos(xx) return 'double MIN(double a, double b) 1100 if a<=b then xx=a else xx=b end if return 'double MAX(double a, double b) 1200 if a>=b then xx=a else xx=b end if return 'double Sign(double a, double b) 1300 if b>=0 then xx=abs(a) else xx=-abs(a) end if return 'double BRENT(double ax, double bx, double cx, double tol, double *xmin) '------------------------------------------------------------------ ' Labels: 1,2 ' 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. '------------------------------------------------------------------ 2000: a=ax : b = cx : gosub 1100 : a=xx 'a=MIN(ax,cx) gosub 1200 : b=xx 'b=MAX(ax,cx) v=bx: w=v: x=v: e=0 xx=x : gosub 1000 : fx=F 'fx=F(x) fv=fx : fw=fx d=0 for iter=1 to ITMAX 'main loop xm=0.5*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.0*tol1 if abs(x-xm)<=tol2-0.5*(b-a) then return 'Test for done here if abs(e)>tol1 then 'Construct a trial parabolic fit r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2#*(q-r) 'bug corrected 07/24/2006 (0.2 instead of 2.0) if q>0 then p=-p q=abs(q) etemp=e e=d if abs(p)>=abs(0.5*q*etemp) or p<=q*(a-x) or p>=q*(b-x) then goto 1 ' 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 then e=a-x else e=b-x end if d=CGOLD*e 2 if abs(d)>=tol1 then u=x+d else a=tol1:b=d:gosub 1300: u=x+xx 'u=x+Sign(tol1,d) end if xx=u : gosub 1000 : fu=F 'fu=F(u) ' This the one function evaluation per iteration if fu<=fx then if u>=x then a=x else b=x end if v=w: fv=fw w=x: fw=fx x=u: fx=fu else if u