' ******************************************************************
' Numerical Gradient
' ******************************************************************

SUB NumGrad(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_ 
            X() AS DOUBLE, Eta AS DOUBLE, G() AS DOUBLE)
' ------------------------------------------------------------------
' Compute the gradient G of function Func at point X using finite
' difference approximations.
' Eta is the relative increment (should not be less than the square 
' root of the machine epsilon ~ about 1.5D-8 in double precision)
' ------------------------------------------------------------------

  CONST EtaMin = 1.5D-8 ' Min. relative increment

  DIM AS INTEGER L, N   ' Array bounds
  DIM AS INTEGER I      ' Loop variable
  DIM AS DOUBLE  Temp   ' Temporary variable
  DIM AS DOUBLE  Delta  ' Increment
  DIM AS DOUBLE  Xm     ' X - Delta
  DIM AS DOUBLE  Xp     ' X + Delta
  DIM AS DOUBLE  Fm     ' F(X - Delta)
  DIM AS DOUBLE  Fp     ' F(X + Delta)

  IF Eta < EtaMin THEN Eta = EtaMin

  L = LBOUND(X)
  N = UBOUND(X)
 
  FOR I = L TO N
    IF X(I) <> 0 THEN
      Delta = Eta * ABS(X(I))
    ELSE
      Delta = Eta
    END IF
 
    Xp = X(I) + Delta
    Xm = X(I) - Delta
 
    Temp = X(I)
    
    X(I) = Xm
    Fm = Func(X())
    
    X(I) = Xp
    Fp = Func(X())
    
    G(I) = (Fp - Fm) / (2 * Delta)
    
    X(I) = Temp
  NEXT I
END SUB

