Adaptive Quadrature: MAT360

Andy Long
Spring, 2024

Here is my version of the adaptive quadrature algorithm from our text, Tea Time Numerical Analysis (Leon Q. Brin): 
	0. Inputs are f, a, b, epsilon, and the maximum depth of recursion. I also input f[a] and f[b]: as we recurse, we can save some computations if we pass values already computed. 
	1. Note that I append the points for each interval to a globally defined list AQpoints (so that I may illustrate the method, or plot the points). That needs to be reset to an empty list before each run.
	2. We use a divide-and-conquer approach: if the error is unacceptable on an interval, then we split the interval in two, split the error tolerance in two, and recurse. 
	3. Maxits keeps track of the depth of the recursion, and so prevents an infinite loop. 
	4. I use the same measure of error as our author.

aTraptiveQuad[f_,enda_,endb_,fenda_,fendb_,eps_,maxits_]:=
Module[{a=enda,b=endb,fa=fenda,fb=fendb,h,c,fc,err,int},
h=(b-a)/2;
c=a+h;
fc=N[f[c]]; (* 
(T2-T1)/3 
= 
h((fa+2fc +fb)/2-(fa+fb))/3 
= 
-h/6(fa-2fc+fb)*)
err=h /6 *(fa -2 fc+ fb); 
If[Abs[err]<eps || maxits == 0
,int=h(fa+4fc +fb)/3;(* Simpson's! Ironic, no? *)
AQpoints=Append[AQpoints,{a,c,b}];
,
int=
aTraptiveQuad[f,a,c,fa,fc,eps/2,maxits-1]
+
aTraptiveQuad[f,c,b,fc,fb,eps/2,maxits-1]
];
Return[int]
]

First I set up the example we’ll use:

f[x_]:=Sin[1/(x+.001)]
a=0;
b=.2;
eps=.001;
its=30;
AQpoints={};

Now we’ll compare Mathematica’s answer and our adaptive algorithm:

NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,N[f[a]],N[f[b]],eps,its]

I stored the points at which we computed the integral, and we’ll see how we can use those same points to make a nice plot:

AQpoints=DeleteDuplicates[AQpoints];
AQpoints=Flatten[AQpoints];
Length[AQpoints]
fvals=f[AQpoints];

Once again I’ll compare Mathematica’s routine and ours:

Plot[f[x],{x,a,b}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,Joined->True
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]

Now we’ll consider the example (sqrt) that our authors provide, p. 298.

f[x_]:=Sqrt[x]
a=0;
b=1;
eps=.00015;
its=30;
AQpoints={};
NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,N[f[a]],N[f[b]],eps,its]
AQpoints=DeleteDuplicates[AQpoints];
AQpoints=Flatten[AQpoints];
Length[AQpoints]
fvals=f[AQpoints];

Plot[f[x],{x,a,b}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,Joined->True
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]

Now we’ll consider the other tough example that the authors provide, p. 299:

f[x_]:=Exp[-100x]
a=0;
b=1;
eps=.00015;
its=30;
AQpoints={};
NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,N[f[a]],N[f[b]],eps,its]
AQpoints=DeleteDuplicates[AQpoints];
AQpoints=Flatten[AQpoints];
Length[AQpoints]
fvals=f[AQpoints];

Plot[f[x],{x,a,b},PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,Joined->True
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]

a=0;
b=0.1;
Show[
Plot[f[x],{x,a,b},PlotRange->{{a,b},{Min[fvals],Max[fvals]}},PlotStyle->{Black,Thick}],
ListPlot[Table[{AQpoints[[i]],fvals[[i]]},{i,Length[AQpoints]}]
,Joined->True,PlotStyle->Yellow
,PlotRange->{{a,b},{Min[fvals],Max[fvals]}}]
]