←to practical programming

Homework "Recursive Adaptive Integration"

(6 points) Recursive open 4-point 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, acc, or relative, eps, accuracy goals.
The integrator has to use a higher order quadrature to estimate the integral and an imbedded (that is, using the same function evaluations) lower order quadrature to estimate the local error.
Something like
```
static double integrate(Func<double,double> f, double a, double b,
double acc=0.001, double eps=0.001, double f2=NaN, double f3=NaN) // NaN indicates first call
{
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 ≤ acc+eps*|Q|) return Q;
else return integrate(f,a,(a+b)/2,acc/√2,eps,f1,f2)+
            integrate(f,(a+b)/2,b,acc/√2,eps,f3,f4);
}
```
Reuse of points is of utmost importance for the effectiveness of the algorithm.
Test your implementation on some interesting integrals, for example (check the expressions before using),
```
∫_₀^¹ dx √(x) = 2/3 ,
∫_₀^¹ dx 1/√(x) = 2 ,
∫_₀^¹ dx √(1-x²) = π/2 ,
∫_₀^¹ 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) = | ²/_√π ∫_₀^{^z} dx exp(-x²)                 , if 0≤z≤1
         | 1-²/_√π ∫_₀^¹ dt exp(-(z+(1-t)/t)²)/t/t , if 1<z
```
make a plot and compare with the tabulated values.
According to a chatbot
```
 erf(1) = 0.84270079294971486934 
```
Now, calculate erf(1) with your routine using eps=0 and decreasing acc=0.1, 0.01, 0.001, … (or something like this). Plot (in log-log scale) the (absolute value of the) difference between your result and the exact result as function of acc.
(3 points) Variable transformation quadratures
• Inplement an (open quandrature) adaptive integrator with the Clenshaw–Curtis variable transformation,
```
 ∫_{_-1}^¹ f(x)dx = ∫_₀^{^π} f(cos(θ))sinθdθ  
```
```
 ∫_{_a}^{^b}dx f(x) = ∫_₀^{^π}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,
```
 ∫_₀^¹ dx 1/√(x) = 2 ,
 ∫_₀^¹ dx ln(x)/√(x) = -4 .
```
• Compare the number of integrand evaluations with the python/numpy's integration routines.
• 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 Clenshaw-Curtis 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.
Hint:
```
int ncalls=0;
Func<double,double> f = z => {ncalls++;return z*z;};
```
(1 point) Error estimate
Make your integrator estimate and return the integration error.
Something like
```
err=|Q-q|
if(err ≤ tol) return (Q,err);
else{
	(Q1,err1)=adapt(...);
	(Q2,err2)=adapt(...);
	return (Q1+Q2,Sqrt(err1*err1+err2*err2));
	}
```
Investigate the quality of this error estimate by calculating some difficult intergrals and comparing the estimated error with the actual error.