src/trapezoidal_and_simpsons_rule_integration/simpson.c

   1 /* integration by simpson rule */
   2 /* M. Vallieres, Spring 2001 */
   3
   4 #include <math.h>
   5
   6 /* function to integrate  */
   7 double fun( double x )
   8 {
   9   double x1, x2, x3, func;
  10
  11   x1 = (x-0.3)/0.1;
  12   x2 = (x-0.5)/0.2;
  13   x3 = (x-0.7)/0.1;
  14   func = exp(-x1*x1)-1.5*exp(-x2*x2)+exp(-x3*x3);
  15   return func;
  16 }
  17
  18 /* simpson rule implementation */
  19 double simpson_int ( double a, double b, int N, double (*f)(double) )
  20 {
  21   double x, dx, integral;
  22   int ix, factor;
  23
  24   if ( N % 2 == 0 )
  25     {
  26       printf( "\n N  must be odd for the Simpson rule \n\n" );
  27       exit(1);
  28     }
  29   dx = (b-a)/(double)(N-1);
  30   integral = f(a) + f(b);
  31   factor = 4;
  32   for ( ix=1 ; ix<N-1 ; ix++ )
  33     {
  34       x = a + dx*ix;
  35       integral = integral + (double)factor*f(x);
  36       factor = 6 - factor;
  37     }
  38   integral = integral * dx / 3.0;
  39   return integral;
  40 }
  41
  42 /* main program */
  43 int main( )
  44 {
  45   double a, b, area;
  46   int N;
  47                           /* integration parameters */
  48   a = 0.0;                /* integration interval */
  49   b = 1.0;
  50
  51   while ( N != 0 )
  52     {
  53                           /* # of grid points */
  54       printf( "\n Input # grid points (0 to end): " );
  55       scanf ( "%d", &N );
  56       if ( N <= 0 )
  57         {
  58           printf( "\n" );
  59           exit(0);
  60         }
  61                           /* calculate integral */
  62       area = simpson_int ( a, b, N, fun );
  63                            /* output results */
  64       printf( " Integral = %15.12f  (Simpson using %d grid points) \n", area, N );
  65
  66     }
  67 }
  68