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
   

   Check that your integrator returns results within the given accuracy goals.

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

            | -erf(-z)                            , if z<0
   erf(z) = | 2/√π ∫0z dx exp(-x²)                 , if 0≤z≤1
            | 1-2/√π ∫01 dt exp(-(z+(1-t)/t)²)/t/t , if 1<z
   

   make a plot, compare with exact results. Can you obtain better accuracy than our approximation from the plot exercise?

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.

   Hint:

   int ncalls=0;
   Func<double,double> f = z => {ncalls++;return z*z;};
   

3. (1 point) Error estimate and 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.