function ShowDefIntegral()
% Illustrate adaptive trapezoidal rule on a function that is very
% nonlinear near zero but which looks more like a horizonal line as
% x gets big.
clc
L = 0.01;
R = 10;
x = linspace(L,R,50000);
semilogx(x,sin(1./x))
shg
tol = .001;
fprintf('Integrate f(x) = sin(1/x) from L = %5.3f to R = %5.3f\n',L,R)
fprintf('Adaptive trapezoidal rule with tol = %10.6f\n\n',tol)
fprintf(' L R R-L\n---------------------------------------\n')
[A,fEvals] = Area1(L,R,@MyF,tol);
fprintf('\n\nIntegral via Area1 = %15.12f Number of f-evals = %4d\n',A,fEvals)
[A,fEvals] = Area2(L,R,@MyF,MyF(L),MyF(R),tol);
fprintf('Integral via Area2 = %15.12f Number of f-evals = %4d\n\n',A,fEvals)

fprintf('\n\n Integral fEvals tol\n------------------------------------\n')
for tol = [.0001 .00001 .000001 .0000001 .00000001]
[A,fEvals] = Area2(L,R,@MyF,MyF(L),MyF(R),tol);
fprintf('%15.12f %5d %10.8f\n',A,fEvals,tol)
end

function [A,fEvals] = Area1(L,R,f,tol)
% Adaptive trapezoidal rule with redundant f-evals
% f is a function of a single real variable that is defined
% on the interval [L,R]
% A is an approximation of the integral of f from L to R.
% fEvals is the number of times that f is called.

% Trapezoidal estimate
h = (R-L);
A1 = h*(f(L)+f(R))/2;
% Trapezoidal estimate based on left and right half intervals...
M = (L+R)/2;
A2 = h*(f(L) + 2*f(M) + f(R))/4;
% Compare and decide whether or not to subdivide the problem...
if abs(A1-A2)<=tol || R-L<=.00001
% Good enough or the interval is so short that we "give up"
A = A2;
fEvals = 5;
fprintf('%10.6f %10.6f %10.6f\n',L,R,R-L)
else
% Recur applying the method to the left and right subproblems
[LeftArea,fEvalsLeft] = Area1(L,M,f,tol);
[RightArea,fEvalsRight] = Area1(M,R,f,tol);
A = LeftArea + RightArea;
fEvals = 5 + fEvalsLeft + fEvalsRight;
end

function [A,fEvals] = Area2(L,R,f,fL,fR,tol)
% Adaptive trapezoidal rule with minimum f-evals
% f is a function of a single real variable that is defined
% on the interval [L,R].
% fL = f(L) and fR = f(R)
% A is an approximation of the integral of f from L to R.
% fEvals is the number of times that f is called.

% Trapezoidal estimate
h = (R-L);
A1 = h*(fL+fR)/2;
% Trapezoidal estimate based on left and right half intervals...
M = (L+R)/2;
fM = f(M);
A2 = h*(fL + 2*fM + fR)/4;
% Compare and decide whether or not to subdivide the problem...
if abs(A1-A2)<=tol || R-L<=.00001
% Good enough or the interval is so short that we "give up"
A = A2;
fEvals = 1;
else
% Recur applying the method to the left and right subproblems...
[LeftArea,fEvalsLeft] = Area2(L,M,f,fL,fM,tol);
[RightArea,fEvalsRight] = Area2(M,R,f,fM,fR,tol);
A = LeftArea + RightArea;
fEvals = 1 + fEvalsLeft + fEvalsRight;
end

function y = MyF(x)
y = sin(1/x);