adaptive_quadrature.c

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

double f(double x, int* nbeval) {
  (*nbeval)++;
  return 1/(1+x*x);
  // replace with 1/(1+x) to obtain log(2) !
}

double intg(double a,double b, double eps, int* nbeval) {
  double x = (b-a)*(f(a, nbeval)+f(b, nbeval))/2;
  if(x>eps) {
    return intg(a,(a+b)/2, eps, nbeval)+intg((a+b)/2,b, eps, nbeval);
  } else {
    return x;
  }
}

int main(void) {
  double eps;
  int nbeval = 0;
  printf("eps\tError\tEvaluations\n");
  for(eps=1; eps > 1e-6; eps /= 10) {
    printf("%f\t%.12lf\t", eps, fabs(intg(0, 1, eps, &nbeval)-.785398163397));
    printf("%d\n", nbeval);
    nbeval = 0;
  }
  return 0;
}