'************************************************************************
'*     Test program: Numerical differentation according to  Romberg     *
'* -------------------------------------------------------------------- *
'* SAMPLE RUN:                                                          *
'* (Find the first derivative of function f(x)=1/x, for x=0.12).        *
'*                                                                      *
'*  Numerical differentation according to Romberg                       *
'*                                                                      *
'*  Test function   f(x) = 1/x                                          *
'*                                                                      *
'*  Put in x-value at which you want to evaluate derivative: 0.12       *
'*  Put in desired accuracy: 1e-10                                      *
'*  Maximal number of columns in Romberg scheme: 4                      *
'*  Starting step size: 0.005                                           *
'*                                                                      *
'*  Results:                                                            *
'*                                                                      *
'*  er_app=  8.893152880773414D-011                                     *
'*  res   = -69.44444444444329                                          *
'*  nend  =  3                                                          *
'*  hend  =  .00062                                                     *
'*                                                                      *
'* -------------------------------------------------------------------- *
'* "Numerical Algorithms with C, By Gisela Engeln-Muellges              *
'*  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                 *
'*                                                                      *
'*                              Basic Release 1.0 By J-P Moreau, Paris. *
'*                                        (www.jpmoreau.fr)             *
'************************************************************************
'Program Test_Difrom

DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

   CLS
   PRINT
   PRINT " Numerical differentation according to Romberg"
   PRINT
   PRINT " Test function   f(x) = 1/x"
   PRINT
   INPUT " Put in x-value at which you want to evaluate derivative: ", x0
   INPUT " Put in desired accuracy: ", eps
   INPUT " Maximal number of columns in Romberg scheme: ", n
   INPUT " Starting step size: ", h

   GOSUB 1000  'call difrom(x,eps,n,h,res,erapp,nend,hend,ierror)

   GOSUB 2000  'call aus(ierror,erapp,res,nend,hend)

END

'test function to derivate
500 'Function func(xx)
  func = 1# / xx
RETURN


1000 'Subroutine difrom(x0, eps, n, h, res, er_app, nend, hend, error)
'*************************************************************************
'*  Computes an approximation for the first derivative of func at x0     *
'*  using the ROMBERG method.                                            *
'* --------------------------------------------------------------------- *
'*  Input parameters :                                                   *
'*                                                                       *
'*    real*8 func(real*8)  (External) Name of function to be             *
'*                         differentiated.                               *
'*    real*8 x0            value of abscissa at which derivative is to   *
'*                         be found.                                     *
'*    real*8 eps           desired accuracy.                             *
'*    integer n            max. number of columns in the Romberg scheme  *
'*                         (n > 1).                                      *
'*    real*8 h             initial step size.                            *
'*                                                                       *
'*  Output parameters:                                                   *
'*                                                                       *
'*    real*8 res           approximate derivative                        *
'*    real*8 er_app        error estimate for  res                       *
'*    integer nend         number of columns actually used in scheme     *
'*    real*8 hend          final step size                               *
'*                                                                       *
'*  Return value (integer ierror):                                       *
'*                                                                       *
'*    0:   no error: er_app < eps                                        *
'*    1:   n < 1  or  eps <= 0  or  h < MACH_EPS                         *
'*    2:   desired accuracy not reached after n steps                    *
'*    3:   step size drooped below  MACH_EPS                             *
'*    4:   lack of sufficient memory (not used here)                     *
'*                                                                       *
'*************************************************************************
  NMAX = 10: XMACHEPS = 1.2D-16
  DIM d(NMAX)

  IF n <= 1 OR eps <= 0# OR h < XMACHEPS THEN
    ierror = 1
    RETURN
  END IF

  h2 = 2# * h
  xx = x0 + h: GOSUB 500: func1 = func
  xx = x0 - h: GOSUB 500
  d(0) = (func1 - func) / h2

'************************************************************************
'* This loop runs until the maximum of Romberg rows is filled or until  *
'* the step size use drops below the machine constant or if the desired *
'* accuracy is reached.                                                 *
'************************************************************************
  ierror = 2
  FOR j = 1 TO n - 1
    d(j) = 0#
    d1 = d(0)
    h2 = h
    h = h * .5#
    nend = j

    IF h < XMACHEPS THEN       'step size less than machine constant
      ierror = 3
      GOTO 10
    END IF

    xx = x0 + h: GOSUB 500: func1 = func
    xx = x0 - h: GOSUB 500
    d(0) = (func1 - func) / h2

    m = 4
    FOR i = 1 TO j
      d2 = d(i)
      d(i) = (m * d(i - 1) - d1) / (m - 1)
      d1 = d2
      m = m * 4
    NEXT i

    erapp = ABS(d(j) - d(j - 1))

    IF erapp < eps THEN
      ierror = 0                         'desired accuracy reached
      GOTO 10
    END IF
  NEXT j

10 res = d(nend)                         'save final values
   hend = h

RETURN


'print results
2000 'Subroutine aus(error, schaetz, res, nend, hend)
  IF ierror <> 1 THEN
    PRINT
    PRINT " Results:"
    PRINT
    PRINT " er_app= "; erapp
    PRINT " res   = "; res
    PRINT " nend  = "; nend
    PRINT " hend  = "; hend
    PRINT
  END IF
RETURN

'end of file tdifrom.bas