'******************************************************** '* 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" [8]. * '* ---------------------------------------------------- * '* 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 * '* * '* BASIC Version By J-P Moreau. * '* (www.jpmoreau.fr) * '******************************************************** defdbl a-h,o-z AX=4.0 : BX=5.0 'X1=4 X2=5 gosub 2000 'Call MNBRAK routine cls print print " The three points are:" print print using " X1=####.###### X2=####.###### X3=####.######"; AX,BX,CX print print " Corresponding function values:" print print using " F1=####.###### F2=####.###### F3=####.######"; FA,FB,FC print print end 'Function to be analyzed 1000 FUNC=xx*SIN(xx)-2.0*COS(xx) 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 'SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,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. 2000 GOLD=1.618034 : GLIMIT=100 : TINY=1e-20 '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. xx=AX:gosub 1000:FA=FUNC 'FA=FUNC(AX) xx=BX:gosub 1000:FB=FUNC 'FB=FUNC(BX) IF FB > FA THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM END IF CX=BX+GOLD*(BX-AX) xx=CX:gosub 1000:FC=FUNC 'FC=FUNC(CX) 1 IF FB >= FC THEN R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) a=ABS(Q-R):b=TINY:gosub 1200:xmax=xx 'xmax=MAX(a,b) a=xmax:b=Q-R:gosub 1300 'xx=Sign(a,b) U=BX-((BX-CX)*Q-(BX-AX)*R)/(2*xx) ULIM=BX+GLIMIT*(CX-BX) IF (BX-U)*(U-CX) > 0 THEN xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) IF FU < FC THEN AX=BX FA=FB BX=U FB=FU GOTO 1 ELSE IF FU > FB THEN CX=U FC=FU GOTO 1 END IF END IF U=CX+GOLD*(CX-BX) xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) ELSEIF (CX-U)*(U-ULIM) > 0 THEN xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) IF FU < FC THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) END IF ELSEIF (U-ULIM)*(ULIM-CX)>=0 THEN U=ULIM xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) ELSE U=CX+GOLD*(CX-BX) xx=U:gosub 1000:FU=FUNC 'FU=FUNC(U) END IF AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GOTO 1 END IF RETURN 'End of file mnbrak.bas