% With the improper integral (integrand improper.m), we found poor convergence
% when the domain was bounded by a singularity. We now address that problem.

O = '---- parabola ---- = O'
Nfor = 1e2;
pow = [1,2,3,4,5,6,7];
Nsteps = (10.^pow).*Nfor
Ntrials = numel(pow) 
for i = 1:Ntrials
   t0 = cputime;
   results(i) = midpoint2(@parabola, 0, pi, Nfor, 10^pow(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 

hold off
plot(log10(Nsteps(1:end-1)), log10(abs(errors)), 'o-k');
xlabel('log10 number of steps', 'fontsize', 16);
ylabel('log10 absolute error', 'fontsize', 16);
title('midpoint rule convergence (domain 0:\pi)', 'fontsize', 18);
hold on


O = '---- improper2 ---- = O'
for i = 1:Ntrials
   t0 = cputime;
   results(i) = midpoint2(@improper2, 0, pi, Nfor, 10^pow(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 
best_estimate = results(end)
error = errors(end)
plot(log10(Nsteps(1:end-1)), log10(abs(errors)), '+-b');


O = '---- improper3 (transformed from improper2) ---- = O'
for i = 1:Ntrials
   t0 = cputime;
   results(i) = midpoint2(@improper3, 0, pi, Nfor, 10^pow(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 
best_estimate = results(end) + 4*sqrt(pi)
   % See derivation in notes. This should match the best estimate of the integral for improper2.
error = errors(end)
O = '---- Calculations complete! ---- = O'

plot(log10(Nsteps(1:end-1)), log10(abs(errors)), '^-g');

legend('parabola','improper2','improper3');

print('integration_test6.pdf','-dpdf')

hold off