Adaptive numerical integration
Example for integral
We want to approximate the integral
with an approximation Q such that Q-I<=1e-2. We want to use the composite trapezoid rule with a subdivision which is adapted to the behavior of the function f: large intervals where |f''| is small, small intervals where |f''| is large.
global Nev % use global variable Nev to count function evaluations
exact = (atan(b^2)-atan(a^2))/2;
xlabel('x'); title('function f(x)')
Try adaptint with tol=.1, .01, .01
Note that the actual error is about 0.3 times smaller than the tolerance.
Try to change the line Q=QT to Q=QS in the code below. What happens with the errors?
figure; fplot(f,[a,b]); hold on
Nev = 2; % 2 initial evaluations at a,b
Nev % number of function evaluations
title(sprintf('tol = %g: adaptint uses %g subintervals',tol,Nev-1))
Recursive function adaptint(f,a,b,Tol)
We use the composite trapezoid rule. On each subinterval we estimate the error |QT-I| by |QT-QS| where QS is the Simpson rule value. We subdivide the intervals until the sum of estimated errors is <= Tol.
In order to demonstrate how this works we
- use the global variable Nev to count the number of function evaluations
- plot a short vertical line at the midpoint to show the resulting mesh
function Q = adaptint(f,a,b,Tol)
QT = (b-a)/2*(fa+fb); % Trapezoid rule
c = (a+b)/2; fc = f(c); % evaluate f in midpoint
Nev = Nev + 1; % count function evaluations
plot([c c],[0,2]*1e-2,'r'); % plot new point
QS = (b-a)/6*(fa+4*fc+fb); % Simpson rule
% for small intervals QS is much closer to I than QT, hence QT-I ~ QT-QS
if abs(QT-QS)<=Tol % if estimated error is <= Tol
Q = QT; % accept trapezoid rule value
Q = adaptint(f,a,c,Tol/2) + adaptint(f,c,b,Tol/2); % use algorithm for [a,c] and [c,b]