'****************************************************************
'* Integration of a discrete function F(x) by Simpson's method  *
'* ------------------------------------------------------------ *
'* SAMPLE RUN                                                   *
'*                                                              *
'* (Find integral of digitalized function sin(x) from x=0 to    *
'*  x=1).                                                       *
'*                                                              *
'* Integration of a discrete real function F(x)                 *
'*             by Simpson's Method                              *
'*                                                              *
'* Number of points: 25                                         *
'* Begin x         : 0                                          *
'* End x           : 1                                          *
'*                                                              *
'* Value of integral: 0.45969770                                *
'*                                                              *
'*                                                              *
'*                               Basic version by J-P Moreau.   *
'*                                    (www.jpmoreau.fr)         *
'****************************************************************
'PROGRAM Test_Discreet_Simpson
DEFDBL A-H, O-Z
DEFINT I, N

NSIZE = 100

DIM Xi(NSIZE), Yi(NSIZE)
'double  dx,result,x,x1,x2

  CLS
  PRINT
  PRINT "    Integration of a discrete real function F(x)"
  PRINT "                by Simpson's Method"
  PRINT
  INPUT " Number of points: ", npoints
  INPUT " Begin x         : ", x1
  INPUT " End x           : ", x2
  'build up test discrete curve F(x)=sin(x)
  dx = (x2 - x1) / (npoints - 1)
  X = x1 - dx
  FOR i = 1 TO npoints
    X = X + dx
    Xi(i) = X: Yi(i) = SIN(X)
  NEXT i

  'call integration function
  GOSUB 1000  'call DiscreetIntegral(npoints,Xi,Yi,result)

  'print value of integral
  PRINT
  PRINT
  PRINT " Value of integral: ";
  PRINT USING "##.########"; result
  PRINT
  PRINT

END 'of main program


'***********************************************************
'* This subroutine returns the integral value of a discrete*
'* real function F(x) defined by n points from x1 to xn.   *
'* ------------------------------------------------------- *
'* INPUTS:                                                 *
'*          n      Number of points (xi, yi)               *
'*          Xi     Table of n xi values in ascending order *
'*          Yi     Table of n yi values                    *
'*                                                         *
'* OUTPUT          The subroutine returns the integral of  *
'*                 F(x) from x1 to xn in result.           *
'* ------------------------------------------------------- *
'* Ref.: "Mathematiques en Turbo Pascal by Alain Reverchon *
'*        and Marc Ducamp, Armand Colin Editor, Paris,     *
'*        1990".                                           *
'***********************************************************
1000 'Subroutine DiscreetIntegral(n, X, Y, result)
'real h,xk,xj,r
  n = npoints
  r = 0#
  i = 1
1010 i1 = i + 1: i2 = i + 2: h = Xi(i1) - Xi(i)
    xk = Xi(i2) - Xi(i1): xj = h + xk
    IF h = 0 THEN
      result = 0#
      RETURN
    END IF
    r = r + ((h + h - xk) * Yi(i) + (xj * xj * Yi(i1) + h * (xk + xk - h) * Yi(i2)) / xk) * xj / 6 / h
    i = i2
  IF (i <= n - 2) THEN GOTO 1010
  IF (n MOD 2) = 0 THEN
    r = r + (Xi(n - 1) - Xi(n - 2)) * (Yi(n - 1) + Yi(n - 2)) / 2#
  END IF
  result = r
RETURN

' end of file disinteg.bas