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 goal δ and realtive accuracy goal ε. The subroutine should return the estimate of the integral, Q, with the error smaller than δ+ε|Q|.

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, δ, or relative, ε, accuracy goals.

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

   Something like

   double integrate(double f(double), double a, double b, double δ, double ε)
   {
   double Q = higher_order_rule, q = lower_order_rule, err=|Q-q|;
   if (err < δ+ε|Q|) return Q;
   else return integrate(f,a,(a+b)/2,δ/√2,ε)+
               integrate(f,(a+b)/2,b,δ/√2,ε);
   }
   

   Reuse points!

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

   ∫01 dx √(x) = 2/3 ,
   ∫01 dx 4√(1-x²) = π
   

2. (3 points) Open quadrature with Clenshaw–Curtis variable transformation

   • Inplement an (open quandrature) adaptive integrator with the Clenshaw–Curtis variable transformation,

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

     ∫abdx f(x) = ∫0πdθ f( (a+b)/2+(b-a)/2*Cos(θ) )*Sin(θ)*(b-a)/2
     

   • Record the number of integrand evaluations and check whether there is any improvement on integrals with integrable divergencies at the end-points of the intervals. For example,

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

   • Calculate the integral

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

     with as many significant digits as possible with and without Clenshaw-Curtis transformation. Compare the accuracy and the number of evaluations.
   • Compare with the GSL's integration routines.

3. (1 point) Infinite limits

   • Make your integrator estimate and return the integration error.
   • 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. Hints: man INFINITY, man isinf.
   • Test your implementation on some (converging) infitine limit integrals and note the number of integrand evaluations.
   • Compare with the GSL's integration routines.