<![CDATA[ PROGRAM AREA
************************************************************************
* Program to approximate the integral of a function over the interval *
* [A,B] using the trapezoidal method. This approximation is calculated *
* by the function subprogram DEFINT; the integrand, the interval of *
* integration, and the # of subintervals are passed as arguments to *
* DEFINT. Identifiers used are: *
* A, B : the endpoints of the interval of integration *
* POLY : the integrand *
* NSUBS : the number of subintervals used *
* *
* Input: A, B, and NSUBS *
* Output: Approximation to integral of F on [A, B] *
************************************************************************
REAL A, B, DEFINT, POLY
INTEGER NSUBS
EXTERNAL POLY
PRINT *, 'ENTER THE INTERVAL ENDPOINTS AND THE # OF SUBINTERVALS'
READ *, A, B, NSUBS
PRINT 10, NSUBS, DEFINT(POLY, A, B, NSUBS)
10 FORMAT (1X, 'APPROXIMATE VALUE USING', I4,
+ ' SUBINTERVALS IS', F10.5)
END
** DEFINT **************************************************************
* Function to calculate the trapezoidal approximation of the integral *
* of the function F over the interval [A,B] using N subintervals. *
* Local variables used are: *
* I : counter *
* DELX : the length of the subintervals *
* X : a point of subdivision *
* Y : the value of the function at X *
* Accepts: Function F, endpoints A and B, and number N of subintervals *
* Returns: Approximate value of integral of F over [A, B] *
************************************************************************
FUNCTION DEFINT(F, A, B, N)
INTEGER N, I
REAL DELX, X, DEFINT, F, A, B
* Calculate subinterval length
* and initialize the approximating sum and X
DELX = (B - A) / REAL(N)
X = A
DEFINT = 0.0
* Now calculate the approximating sum
DO 10 I = 1, N - 1
X = X + DELX
Y = F(X)
DEFINT = DEFINT + Y
10 CONTINUE
DEFINT = DELX * ((F(A) + F(B)) / 2.0 + DEFINT)
END
**POLY******************************************************************
* The integrand *
************************************************************************
FUNCTION POLY(X)
REAL X, POLY
POLY = X ** 2 + 1
END
]]>