% script to perform adaptive quadrature
clear all, close all global pts % function to be integrated defined in routine f
f = 'integrand';
a = 1; b = 3;
pts = [a;b]; tol = input('Enter error tolerance: '); % this is just to plot the graph
% it's usually a good idea to look at the integrand
% if possible before you try to integrate it
ezplot(f,[1,3])
disp('Hit any key to continue')
pause
hold on fa = feval(f,a);
fb = feval(f,b); Sf_old = simp(a,b,f,fa,fb); Sf = adaptiveSimpson(a,b,f,fa,fb,tol,Sf_old) qpts = sort(pts);
plot(qpts,zeros(length(pts),1),'rx') disp('number of function evaluations')
disp(length(pts))
disp('Hit any key to continue')
pause % now compare result with straight Simpson's rule
% using the same number of points
sum = 0;
h = (b-a)/(length(pts)-1);
for i=0:length(pts)-1,
fxi = feval(f,a+i*h);
if i == 0 | i == length(pts)-1,
sum = sum + fxi;
elseif mod(i,2) == 1,
sum = sum + 4*fxi;
else
sum = sum + 2*fxi;
end
end
disp('Simpson''s rule with the same number of points')
sum = h/3*sum % compute exact solution
% anti-derivative of integrand is 10*cos(10/x)
% so ...
exactSolution = 10*(cos(10/b)-cos(10/a)); errorAdaptiveSimpson = exactSolution - Sf
errorUniformSimpson = exactSolution - sum