Homework "Adaptive Integration"

1. (6 points) Recursive adaptive integrator

   Implement a recursive open-quadrature 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

   static double integrate
   (Func<double,double> f, double a, double b,
   double δ=0.001, double ε=0.001, double f2=NaN, double f3=NaN)
   {
   double h=b-a;
   if(IsNaN(f2)){ f2=f(a+2*h/6); f3=f(a+4*h/6); } // first call, no points to reuse
   double f1=f(a+h/6), f4=f(a+5*h/6);
   double Q = (2*f1+f2+f3+2*f4)/6*(b-a); // higher order rule
   double q = (  f1+f2+f3+  f4)/4*(b-a); // lower order rule
   double err = |Q-q|;
   if (err ≤ δ+ε|Q|) return Q;
   else return integrate(f,a,(a+b)/2,δ/√2,ε,f1,f2)+
               integrate(f,(a+b)/2,b,δ/√2,ε,f3,f4);
   }
   

   Remember to reuse points!

   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 4√(1-x²) = π ,
   ∫01 dx ln(x)/√(x) = -4
   

   Using your integrator implement the error function via its integral representation,

   erf(z) = ∫0z dx 2/√π exp(-x²) ,
   

   make a plot, compare with exact results.

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
   

   Calculate some integrals with integrable divergencies at the end-points of the intervals; record the number of integrand evaluations; compare with your ordinary integrator without variable transformation. For example,

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

   Compare the number of integrand evaluations with the python/numpy'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 integrator.

   Hints: double.PositiveInfinity, IsInfinity, IsNegativeInfinity, IsPositiveInfinity.

   Test your implementation on some (converging) infitine limit integrals and note the number of integrand evaluations.

   Compare with the python/numpy integration routines.