/* Examples using trapezoidal rule for integration */

#include        <stdio.h>
#include        <math.h>

double trapezoidal(double (*f)(double),
                   double a, double b, int n)
{
  double h, sum;
  int    i;
  h = (b-a)/n;
  sum = 0.5*((*f)(a)+(*f)(b));
  for ( i = 1; i < n; i++ )
    sum = sum + (*f)(a+i*h);
  return sum*h;
}


double x_sqr_cosx(double x)
{  return x*x*cos(M_PI*x);  }


int main(int argc, char *argv[])
{
  int        i, n;
  double h, sum, sum1, sum2, sum3, x;
  double fval1, fval2;

  /* trap1.c */
  n = 100;
  h = 1.0/n;
  sum = (exp(0.0)+exp(1.0))*(h/2);
  for ( i = 1; i < n; i++ )
    {
      x = i*h;
      sum = sum + exp(x)*h;
    }
  
  /* trap2.c */
  n = 200;
  h = 1.0/n;
  sum2 = (exp(0.0)+exp(1.0))*(h/2);
  for ( i = 1; i < n; i++ )
    {
      x = i*h;
      sum2 = sum2 + exp(x)*h;
    }
  printf("Error is approximately %g\n", sum2-sum);
  
  /* trap3.c */
  n = 100;
  h = 1.0/n;
  x = 0.0; fval1 = x*x*cos(M_PI*x);
  x = 1.0; fval2 = x*x*cos(M_PI*x);
  sum = (fval1+fval2)*(h/2);
  for ( i = 1; i < n; i++ )
    {
      x = i*h;
      sum = sum + x*x*cos(M_PI*x)*h;
    }
  
  /* trap4.c */
  sum  = trapezoidal(exp,0.0,1.0,100);
  sum2 = trapezoidal(exp,0.0,1.0,200);
  printf("Error is approximately %g\n", sum2-sum);
  /* trap5.c */
  sum = trapezoidal(x_sqr_cosx,0.0,1.0,100);
  
  printf("integral of x^2.cos(pi*x), x = 0..1 approx= %g\n", sum);
}