%% Adaptive numerical integration
%% Example for integral
% We want to approximate the integral
% I = \int_0^4 \frac{x}{1+x^4} dx
% 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.

clearvars; 
global Nev                    % use global variable Nev to count function evaluations
f = @(x) x./(1+x.^4);
a = 0; b = 4;
exact = (atan(b^2)-atan(a^2))/2;
fplot(f,[a,b])
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?

for tol = [.1 .01 .001]
  figure; fplot(f,[a,b]); hold on
  Nev = 2;                       % 2 initial evaluations at a,b
  Q = adaptint(f,a,b,tol)
  Nev                            % number of function evaluations
  error = Q - exact              
  title(sprintf('tol = %g: adaptint uses %g subintervals',tol,Nev-1))
  hold off;
end

%
%% 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)
  global Nev             
  fa = f(a); fb = f(b);
  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
  else                         
    Q = adaptint(f,a,c,Tol/2) + adaptint(f,c,b,Tol/2);  % use algorithm for [a,c] and [c,b]
  end
end