Numerical Integration

Trapezoid method

/* trapezoid.c */

#include <stdio.h>

double f( double x) {

  return(x*x) ;
}

int main() {

  double a, b, delta, integral;
  int i,n;

  a = 0;
  b = 1;
  n = 100000;

  delta = (b-a)/n;
  integral = 0.0;

  for (i=0; i < n; i=i+1)
    integral = integral + (f(a + i*delta)+f(a+(i+1)*delta))/2*delta;

  printf("The integral of f from %f to %f is equal to %f\n",a,b,integral);
}   

Simpson's method

/* simpson.c */


#include <stdio.h>

double f(double x) {

  return(x*x) ;
}

int main() {

  double a, b, delta, integral;
  int i,n;

  a = 0;
  b = 1;
  n = 100000;

  delta = (b-a)/n;
  integral = 0.0;

  for (i=0; i < n; i=i+1)
    integral = integral + (f(a+i*delta)+4*f(a+(i+0.5)*delta)+f(a+(i+1)*delta))/6*delta;

  printf("The integral of f from %lf to %lf is equal to %lf\n",a,b,integral);
}   

Monte Carlo Method

#include <stdio.h>

#include <stdlib.h>
#define RANDOM_MAX 2147483647

int main() {
 
    int i=0, j=0;
    float x, y;

    while(1) {
      x = (double)random()/(double)RANDOM_MAX;
      y = (double)random()/(double)RANDOM_MAX;
    if ((x*x+y*y) < 1)
       j++;
    i++;
   if (i%1000000 == 1)
   printf("%f\n",4*(float)j/(float)i);
   }
}