!******************************************************** !* 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.d0 X2=5.d0 * !* * !* 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 * !* * !******************************************************** PROGRAM TEST_MNBRAK REAL*8 X1, X2, X3, F1, F2, F3 X1=4.D0; X2=5.D0 CALL MNBRAK(X1,X2,X3,F1,F2,F3) print *,' ' print *,' The three points are:' print *,' ' write(*,10) X1, X2, X3 print *,' ' print *,' Corresponding function values:' print *,' ' write(*,20) F1, F2, F3 print *,' ' print *,' ' stop 10 format(' X1=',F10.6,' X2=',F10.6,' X3=',F10.6) 20 format(' F1=',F10.6,' F2=',F10.6,' F3=',F10.6) end !Function to be analyzed REAL*8 FUNCTION FUNC(X) REAL*8 X FUNC=X*SIN(X)-2.D0*COS(X) RETURN END 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. IMPLICIT REAL*8 A-H,O-Z PARAMETER(GOLD=1.618034,GLIMIT=100.,TINY=1.D-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. FA=FUNC(AX) FB=FUNC(BX) IF(FB.GT.FA) THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF CX=BX+GOLD*(BX-AX) FC=FUNC(CX) 1 IF(FB.GE.FC) THEN R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) IF((BX-U)*(U-CX).GT.0) THEN FU=FUNC(U) IF(FU.LT.FC) THEN AX=BX FA=FB BX=U FB=FU GOTO 1 ELSE IF(FU.GT.FB) THEN CX=U FC=FU GOTO 1 ENDIF U=CX+GOLD*(CX-BX) FU=FUNC(U) ELSE IF((CX-U)*(U-ULIM).GT.0) THEN FU=FUNC(U) IF(FU.LT.FC) THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=FUNC(U) ENDIF ELSE IF((U-ULIM)*(ULIM-CX).GE.0) THEN U=ULIM FU=FUNC(U) ELSE U=CX+GOLD*(CX-BX) FU=FUNC(U) ENDIF AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GOTO 1 ENDIF RETURN END !End of file mnbrak.f90