'********************************************************
'*    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