Exercise Adaptive Integration

Introduction:

In this exercise you have to implement a subroutine which integrates a given function f(x) from a to b with the absolute accuracy acc and realtive accuracy eps. The subroutine should return the estimate of the integral, Q 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;
}

This 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 my_function inside the main function (a technique called 'nested functions') where it can inherit all variables from the main-function scope.

Exercises

1. (6 points) Recursive adaptive integrator

   Implement a recursive adaptive integrator that estimates the integral of a given function f(x) on a given interval [a,b] with the required absolute, acc, or relative, eps, accuracy goals.

   The integrator has to use a combination of a higher order quadrature and also an imbedded lower order quadrature to estimate the local error.

   Reuse points!

   The integrator has to use either closed quadratures (which use the end-points and the middle-point of the interval) or open quadratures (which use neither the end-points nor the middle-point of the interval).

   Test your implementation on some interesting integrals, for example (check the expressions before using),

   ∫01 dx √(x) = 2/3 ,
   ∫01 dx 1/√(x) = 2 ,
   ∫01 dx ln(x)/√(x) = -4 .
   

   Calculate

   ∫01 dx 4√(1-(1-x)2) = π
   

   with as many significant digits as you can and count the number of integrand evaluations.

2. (3 points) Clenshaw–Curtis variable transformation

   Combine your adaptive integrator with the Clenshaw–Curtis variable transformation,

   ∫-11 f(x)dx = ∫0π f(cos(θ))sinθdθ  
   

   Check whether there is any improvements on integrals with integrable divergencies at the ends of the intervals.

   Compare with library routines.

3. (1 point) Infinite limits.

   Generalize your integrator to 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. You might need to use the

   -std=c99

   switch.

   In Python3 import math and you will have math.inf.

   Test your implementation on some (converging) infitine limit integrals and note the number of integrand evaluations. Compare with library routines (for example, gsl_integration_qags).