//#################################################################### // Estimating PI using the fact that the integral of 1/(1+x*x) from 0 // to 1 is PI/4. // The purpose of the program is to compare 3 different numerical // integration methods. //#################################################################### #include typedef double (*DfD) (double); double midpoint_int (DfD f, double x0, double x1, int n) { int i; double x, dx, sum = 0.0; dx = (x1-x0)/ n; for (i = 0, x = x0 + dx/2; i < n; i++, x += dx) sum += f(x); return sum * dx; } double trapezoidal_int (DfD f, double x0, double x1, int n){ double x, dx, sum; int i; dx = (x1-x0)/ n; sum = (f(x0) + f(x1))/2; for (i=1,x = x0 + dx; i < n; i++, x += dx) sum += f(x); return sum * dx; } double simpsons_int(DfD f, double x0, double x1, int n){ double x, sum, dx = (x1 - x0) / n; sum = f(x1) - f(x0); for(x = x0; x+dx/2 < x1; x += dx) sum += 2.0 * f(x) + 4.0 * f(x + dx/2); return sum * dx / 6.0; } double fun(double x) { return 1.0/(1+x*x); } int main() { double x0=0, x1=1; int n=100; double mp = 4*midpoint_int(fun, x0, x1, n); double tz = 4*trapezoidal_int(fun, x0, x1, n); double si = 4*simpsons_int(fun, x0, x1, n); printf("midpoint: pi= %f\n", mp); printf("trapezoidal: pi= %f\n", tz); printf("simpson: pi= %f\n", si); }