Exercise Adaptive Integration

Hints:

In this exercise you will have to implement a subroutine which integrates a given function f from a to b with the absolute accuracy acc and realtive accuracy eps. The subroutine should return the estimate of the integral, Q, and the estimate of the error, err, such that the following conditions are satisfied:

• err ~ |Q - ∫ab f(x)dx| ,
• err < acc + |Q|*eps .

In the C language the interface to the subroutine could be

• double my_integrator(
  double f(double), double a, double b, double acc, double eps, double *err)
  {
  Q = my_estimate_of_the_integral;
  *err = my_estimate_of_the_error; // should, of course, be smaller than acc+|Q|*eps
  return Q;
  }
  

That is, the subroutine returns the estimate of the integral and writes the estimate of the error into the provided variable double *err.
The main function could than call the integrator as

• int main (void)
  {
  double a=0, b=1, err, acc=0.0001, eps=0.0001;
  double my_function(double x){return x*x;}
  double Q=my_integrator(my_function, a, b, acc, eps, &err); // remember!: &err
  printf("integral=%g, error=%g\n",Q,err);
  return 0;
  }
  

In C++ language—in addition to the C interface where the user supplies the pointer to his function and the pointer to his err variable—one can also use the std::function object; and pass the err variable by reference:

• double my_CXX_integrator(
  	std::function<double(double)> f,
  	double a, double b, double acc, double eps, double& err)
  

Exercises

1. (6 points) Recursive adaptive integrator
   • Implement a recursive adaptive integrator which estimates the integral of a given function f on a given interval [a,b] with the required absolute, acc, or relative, eps, accuracy.
   • Test your implementation on some interesting integrals.
   • Calculate ∫01  dx (ln(x)/√(x)) = -4 with acc=eps=0.001 and estimate the number of integrand evaluations. Since the integrand is singular at the left end-point of the integration interval, you must use open quadratures – the quadratures which do not use the end-points of the interval. For example,
xi = { 1/6, 2/6, 4/6, 5/6 } (nodes)
wi = { 2/6, 1/6, 1/6, 2/6 } (trapezium rule)
vi = { 1/4, 1/4, 1/4, 1/4 } (rectangle rule)

   • Calculate ∫01  dx 4√(1-(1-x)2) = π with many significant digits and estimate the number of integrand evaluations. If you have implemented both closed and open quadratures, compare the effectiveness of the two.
2. (3 points) Variable transformations
   • Generalize your integrator such that it can accept infinite limits. An infinite limit integral can be converted by a variable transformation (see lecture notes) into a finite limit integral, which can then be evaluated by your ordinary integrator. In C/C++ infinities are defined by macros INFINITY and ‑INFINITY and can be identified by the macro isinf. Test it on some (converging) infitine limit integrals.
   • Make a variant of your integrator which implements the Clenshaw-Curtis variable transformation method,
     ∫-11f(x)dx = ∫0π f(cosθ) sinθ dθ ,
which should cope with certain integrable singularities. Test it on some interesting integrable singularities.
     Hint:
     ∫abf(x)dx = ∫-11 f((a+b)/2 + (b-a)/2 t) (b-a)/2 dt ,

3. (1 point) Comparison with ODE driver
   • A definite integral Q=∫abf(x)dx can be reformulated as an ordinary differential equation, y'=f(x), y(a)=0, y(b)=Q, which can be solved with your adaptive ODE solver. Pick an interesing f(x) and compare the effectiveness of your ODE drivers with your adaptive integrators.
4. (0 points)
   • Devise some higher order classical quadratures and implement them in your adaptive integrator.
   • Try double Clenshaw-Curtis transformation and see whether it can even better cope with integrable singularities