Adaptive Quadrature: MAT360

Andy Long
Spring, 2015

Here is my version of the adaptive quadrature algorithm from our text:
  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. This should really be an optional argument, but I don’t know how to do that with Mathematica!
  1. Note that I append the points for each interval to a globally defined list AQpoints (so that I may illustrate the method). 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 authors.

adaptiveQuadrature[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]];
err=h /6 *(fa -2 fc+ fb);
If[Abs[err]<eps || maxits == 0
,int=3*err+2*h*fc;AQpoints=Append[AQpoints,{a,c,b}];
,int=adaptiveQuadrature[f,a,c,fa,fc,eps/2,maxits-1]
+
adaptiveQuadrature[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}]
adaptiveQuadrature[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}]
adaptiveQuadrature[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 our authors provide, p. 299:

f[x_]:=Exp[-100x]
a=0;
b=1;
eps=.00015;
its=30;
AQpoints={};
NIntegrate[f[x],{x,a,b}]
adaptiveQuadrature[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]}}]