/* % % Written by: % -- % John L. Weatherwax 2005-05-26 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- Recursive trapezoidal rule (serial). */ #include #include #include /* Global storage for the integral: */ double sum = 0.0; /* Count the number of times we recurse: */ int rec_indx = 0; /* The function to integrate: */ double f( double x ){ return 4/(1+x*x); /* return sin(x)*exp(x); */ /* return 1.0; */ } double rec_quad( double tol, double xL, double fL, double xR, double fR, double area0 ){ double xM,fM,aL,aR; /* 1) Split the interval into two 2) compute two areas with the trapezoidal rule 3) Stop if the relative tolerance is met */ rec_indx = rec_indx+1; xM = 0.5*(xL+xR); fM = f(xM); aL = 0.5*(xM-xL)*(fM+fL); aR = 0.5*(xR-xM)*(fM+fR); printf( "rec_quad #(%d): aL=%f, aR=%f, aL+aR=%f, area0=%f\n", rec_indx, aL, aR, aL+aR, area0 ); printf( " : fabs(area0-(aL+aR))=%f; fabs(area0-(aL+aR))/area0=%f\n", fabs(area0-(aL+aR)), fabs(area0-(aL+aR))/area0 ); if( fabs(area0-(aL+aR))/area0 <= tol ) return aL+aR; else return rec_quad( tol, xL, fL, xM, fM, aL ) + rec_quad( tol, xM, fM, xR, fR, aR ); } int main( int argc, char* argv[] ){ double tol,xL,fL,xR,fR,area0; /* Argument checking */ if( argc != 4 ){ fprintf( stderr, "Expecting three input arguments\n" ); exit(1); } /* Input the left/right end points of the inteval */ xL = atof(argv[1]); xR = atof(argv[2]); /* Input the tolerance required: */ tol = atof(argv[3]); fL = f(xL); fR = f(xR); area0 = 0.5*(xR-xL)*(fL+fR); printf( "main: fL=%f, fR=%f, area0=%f\n", fL, fR, area0 ); sum = rec_quad( tol, xL, fL, xR, fR, area0 ); /* Display the resulting integral: */ printf("Integral = %18.10f\n", sum); }