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 * * SIN : the integrand (library function) * * NSUBS : the number of subintervals used * * * * Input: A, B, and NSUBS * * Output: Approximation to integral of F on [A, B] * ************************************************************************ REAL A, B, DEFINT INTEGER NSUBS INTRINSIC SIN PRINT *, 'ENTER THE INTERVAL ENDPOINTS AND THE # OF SUBINTERVALS' READ *, A, B, NSUBS PRINT 10, NSUBS, DEFINT(SIN, 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