Exercise Adaptive Integration

Introduction:

In this exercise you 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 condition is satisfied:

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

The subroutine might also return the estimate of the error

err ≈ | Q - ∫ab f(x)dx | 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 be smaller than acc+|Q|*eps */
  	return Q;
  }
  

where the subroutine returns the estimate of the integral and stores the estimate of the error in the provided variable double *err.
The main function could then call the integrator as

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

The gcc compiler actually allows the definition of the my_function inside the main function (a technique called 'nested functions') where it can inherit all variables from the main-function scope.

In C++ language — in addition to the C-interface where the user supplies the pointer to their function and the pointer to their 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)
  

In Java the integrator usually defines a public interface,

• public static interface integrand{ public double call(double x); }
  

and then the user creates an object which implements this interface and passes the object as a parameter to the integrator.

Exercises

1. (6 points) Recursive adaptive integrators

   • Implement two recursive adaptive integrators which estimate the integral of a given function f on a given interval [a,b] with the required absolute, acc, or relative, eps, accuracy goals.
     • One integrator has to use a closed quadratures (which use the end-points of the interval) – a combination of a higher order closed quadrature and an imbedded lower order closed quadrature (for error estimation). Remember to reuse points.
     • The other integrator has to use open quadratures (which do not use the end-points of the interval) – a combination of a higher order open quadrature and an imbedded lower order open quadrature (for error estimation). Reuse points.
   • Test your implementations on some interesting integrals and compare the number of integrand calls for your integrators.
   • Calculate ∫01  dx 4√(1-(1-x)2) = π with as many significant digits as you can and estimate the number of integrand evaluations. Compare with library routines (for example, gsl_integration_qags).

2. (3 points)

   • Now re-implement your adaptive integrator without the recursive call to itself.

   • Generalize your integrators such that they 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 in math.h by macros INFINITY and ‑INFINITY and can be identified by the macro isinf. Test your implementation on some (converging) infitine limit integrals.

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, Q=y(b), 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 point) Singular integrands

   • Devise a quadrature which can integrate exactly 1/√x and/or log(x) etc.
   • Implement an adaptive integrator which generally uses closed quadrature but switches to this quadrature in the vicinity of a detected singularity.
   • Calculate ∫01  dx (ln(x)/√(x)) = -4 with acc=eps=0.001 and estimate the number of integrand evaluations with you open-adaptive integrator from exercise A and the present integrator.
   • Compare with library routines (for example, gsl_integration_qags).