Hallo!
This code computes a double integral using the
Trapezoidal Rule:
PHP Code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Evaluate the required double integral by the use of the Trapezoidal Rule: */
double trap ( double aa, /* Lower bound of x's interval. */
double bb, /* Higher bound of x's interval. */
double cc, /* Lower bound of y's interval. */
double dd, /* Higher bound of y's interval. */
int NX, /* Number of point to discretize x's interval. */
int NY); /* Number of point to discretize y's interval. */
/* Evaluate the required integrand of the double integral, which at the same
time is another integral that is going to be evaluated with the Trapezoidal
Rule: */
double feval ( double xx, /* Function argument (mathematically speaking). */
double hh, /* Step size for local Trapezoidal Rule. */
double cc, /* Lower bound of y's interval. */
double dd, /* Higher bound of y's interval. */
int NN); /* Number of point to discretize y's interval. */
/* Evaluate the actual function F(xx, yy) that stands for the
"integrand of the integrand" of the double integral. This is also useful
because we can easily try for different test cases just by changing this
module. */
double IntegrandFF (double xx, double yy);
int main (void) {
double aa; /* Lower bound of x's interval. */
double bb; /* Higher bound of x's interval. */
double cc; /* Lower bound of y's interval. */
double dd; /* Higher bound of y's interval. */
double result; /* Result of the entire integration. */
int NX; /* Number of point to discretize x's interval. */
int NY; /* Number of point to discretize y's interval. */
/* Parameters of the problem: */
aa = 0.0;
bb = M_PI;
cc = 0.0;
dd = M_PI;
NX = 10000;
NY = 10000;
result = trap (aa, bb, cc, dd, NX, NY);
printf("The resulting double integration equals %lf\n", result);
return EXIT_SUCCESS;
}
/* Evaluate the actual function F(xx, yy) that stands for the
"integrand of the integrand" of the double integral. This is also useful
because we can easily try for different test cases just by changing this
module. */
double IntegrandFF(double xx, double yy) {
return sin(xx)*sin(yy);
}
/* Evaluate the required integrand of the double integral, which at the same
time is another integral that is going to be evaluated with the Trapezoidal
Rule: */
double feval (double xx, double hy, double cc, double dd, int NN) {
double yy; /* Spatial variable per each iteration. */
double acum_01d; /* Auxiliary accumulation variable. */
unsigned int ii; /* Iterator. */
yy = cc;
acum_01d = 0.0;
for (ii = 1; ii < NN; ii++) {
yy += hy;
acum_01d += IntegrandFF(xx, yy);
}
return hy*(0.5*(IntegrandFF(xx, cc) + IntegrandFF(xx, dd)) + acum_01d);
}
/* Evaluate the required double integral by the use of the Trapezoidal Rule: */
double trap (double aa, double bb, double cc, double dd, int NX, int NY) {
double hx; /* Grid spacing for x. */
double hy; /* Grid spacing for y. */
double xx; /* Spatial variable per each iteration. */
double acum_01d; /* Auxiliary accumulation variable. */
unsigned int ii; /* Iterator. */
printf("We will integrate from %f to %f and from %f to %f\n", aa, bb, cc, dd);
printf("We will integrate using %i points for x and %i for y\n", NX, NY);
/* Compute the grid spacing values: */
hx = (bb - aa)/NX;
hy = (dd - cc)/NY;
printf("That yields to hx = %f and hy = %f\n", hx, hy);
xx = aa;
acum_01d = 0.0;
for (ii = 1; ii < NX; ii++) {
xx += hx;
acum_01d += feval(xx, hy, cc, dd, NY);
}
return hx*(0.5*(feval(aa, hy, cc, dd, NY) + feval(bb, hy, cc, dd, NY))
+ acum_01d);
}
I would like to time the execution of the code using
time.h.
Specifically the time the function
trap() takes to execute. What would be better? To start/finish the timing outside the routine or to do it inside of it?
I would like to time it from the source code, i.e. not using external software such as the UNIX function
time or profilers.
Also, I would like to time it in seconds with precission of milliseconds.
Any other related comments to timing code in ANSI C?
Please, feel free to "destroy" my code, in the sense of criticize it!
Best regards to everyone!
Timing a code: what's better?
