'*********************************************************
'*  Estimate the nth derivative of a real function f(x)  *
'*  at point x with step h (n from 1 to 5).              *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* Nth derivative of a real function f(x):               *
'*                                                       *
'* Function to derivate: x^3                             *
'*                                                       *
'* At point x0= 10                                       *
'* Order of derivative (1 to 5): 1                       *
'* Value of step h= 0.1                                  *
'*                                                       *
'* f"(x0)=300.000000                                     * 
'*                                                       *
'* At point x0= 10                                       *
'* Order of derivative (1 to 5): 2                       *
'* Value of step h= 0.1                                  *
'*                                                       *
'* f"(x0)= 60.000000                                     *
'*                                                       *
'* ----------------------------------------------------- *
'* Ref.: "Analyse numerique en C By Alain Reverchon and  *
'*        Marc Ducamp, Armand Colin Editeur, Paris,      *
'*        1990" [BIBLI 14].                              *
'*                                                       *
'*                          Basic Version By J-P Moreau. *
'*                               (www.jpmoreau.fr)       *
'*********************************************************
'Program Nderiv

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

'Main variables:
'double  d: nth derivative at x
'        h: increment step
'        x: abscissa of function f(x)
'integer n: order of derivative (1 to 5)

  CLS
  PRINT
  PRINT " Nth derivative of a real function f(x):"
  PRINT
  PRINT " Function to derivate: x^3"
  PRINT
  INPUT " At point x0= ", x
  INPUT " Order of derivative (1 to 5): ", n
  INPUT " Value of step h= ", h
  PRINT

  GOSUB 1000  'd = nderiv(x,n,h)

  IF n = 1 THEN
    PRINT USING " f'(x0) = ####.######"; d
  ELSEIF n = 2 THEN
    PRINT USING " f''(x0) = ####.######"; d
  ELSEIF n = 3 THEN
    PRINT USING " f'''(x0) = ####.######"; d
  ELSEIF n = 4 THEN
    PRINT USING " f''''(x0) = ####.######"; d
  ELSEIF n = 5 THEN
    PRINT USING " f'''''(x0) = ####.######"; d
  END IF

  IF n < 1 OR n > 5 THEN
    PRINT " n must be between 1 and 5 !"
  END IF

  PRINT

END

'sample function f(xx)=xx^3
500 f = xx * xx * xx
RETURN

'---------------------------------------------------
' returns the value of nth derivative of function
' f(x) at point x0 with increment step h
' Note that for polynomial functions h must not be
' too small. Here h=1 gives same results than h=0.01
'---------------------------------------------------
1000 'Function nderiv(x, n, h)

DIM itable(5, 7)

  itable(1, 1) = 3: itable(1, 2) = 45: itable(1, 3) = -9: itable(1, 4) = 1
  itable(1, 5) = 0: itable(1, 6) = 0: itable(1, 7) = 60
  itable(2, 1) = 3: itable(2, 2) = 270: itable(2, 3) = -27: itable(2, 4) = 2
  itable(2, 5) = 0: itable(2, 6) = 0: itable(2, 7) = 180
  itable(3, 1) = 4: itable(3, 2) = -488: itable(3, 3) = 338: itable(3, 4) = -72
  itable(3, 5) = 7: itable(3, 6) = 0: itable(3, 7) = 240
  itable(4, 1) = 4: itable(4, 2) = -1952: itable(4, 3) = 676: itable(4, 4) = -96
  itable(4, 5) = 7: itable(4, 6) = 0: itable(4, 7) = 240
  itable(5, 1) = 5: itable(5, 2) = 1938: itable(5, 3) = -1872: itable(5, 4) = 783
  itable(5, 5) = -152: itable(5, 6) = 13: itable(5, 7) = 288

  IF n > 5 THEN
    d = 0#
    RETURN
  END IF

  xx = x: GOSUB 500: fa = f   'fa=f(x)
  d = 0#
  FOR i = 1 TO itable(n, 1)
    xx = x - i * h: GOSUB 500: fb = f    'fb=f(x-i*h)
    xx = x + i * h: GOSUB 500: fc = f    'fc=f(x+i*h)
    IF (n MOD 2) <> 0 THEN
      d = d + itable(n, i + 1) * (fc - fb)
    ELSE
      d = d + itable(n, i + 1) * ((fc - fa) + (fb - fa))
    END IF
  NEXT i
  d = d / itable(n, 7)
  FOR i = 1 TO n
    d = d / h
  NEXT i

RETURN

'end of file nderiv.bas