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