Computational Physics Lectures: Numerical integration, from Newton-Cotes quadrature to Gaussian quadrature

Code example

A simple function which implements this algorithm is as follows

double TrapezoidalRule(double a, double b, int n, double (*func)(double))
{
      double TrapezSum;
      double fa, fb, x, step;
      int    j;
      step=(b-a)/((double) n);
      fa=(*func)(a)/2. ;
      fb=(*func)(b)/2. ;
      TrapezSum=0.;
      for (j=1; j <= n-1; j++){
         x=j*step+a;
         TrapezSum+=(*func)(x);
      }
      TrapezSum=(TrapezSum+fb+fa)*step;
      return TrapezSum;
}  // end TrapezoidalRule

The function returns a new value for the specific integral through the variable TrapezSum.