% ### EXintegration1.m ###  2016.01.13 CB

% Numerical integration example 

% (modified from) Original source:
% http://ef.engr.utk.edu/ef230-2011-01/modules/matlab-integration/

clear;
% ----------------------
% User parameters
F = @(x)(sin(x)); % function to integrate
%F = @(x)(exp(-x.^2/2)); % function to integrate
xL= [0 pi]; % integration limits

N= 2;    % Method A - # of "rectangles" for LEFT and RIGHT
pts= [3 4 5 6 10 25]; % Method B - # of points to consider integrating (via trapz function)
dur= 1;     % Method B - pause duration [s] for trapz loop
% ----------------------

% ***************
% Show the curve 
figure(1);
fplot(F,[xL(1),xL(2)]) % a quick way to plot a function
xlabel('x'); ylabel('F(x)');

% ***************
% Method A
% Approximate the integral via brute force LEFT and RIGHT Riemann sums
sumL= 0; sumR=0; 
delX= (xL(2)-xL(1))/N;    % step-size
x= linspace(xL(1),xL(2),N+1);  % add one since N is # of 'boxes' and is really N-1
for nn=1:N
    sumL= sumL + F(x(nn))*delX;
    sumR= sumR + F(x(nn+1))*delX;
end
disp(['left-hand rule yields =',num2str(sumL),' (for ',num2str(N+1),' steps)']);
disp(sprintf('right-hand rule yields = %g', sumR));
disp(['area calculated by 1/2*(LEFT+RIGHT) for ',num2str(N),' points =',num2str(0.5*(sumL+sumR))]);

% ***************
% Method B
% Approximate the integral via trapz for different numbers of points
for np=pts
    figure(2); clf % clear the current figure
    hold on % allow stuff to be added to this plot
    x = linspace(xL(1),xL(2),np); % generate x values
    y = F(x); % generate y values
    a2 = trapz(x,y); % use trapz to integrate
    % Generate and display the trapezoids used by trapz
    for ii=1:length(x)-1
        px=[x(ii) x(ii+1) x(ii+1) x(ii)];   py=[0 0 y(ii+1) y(ii)];
        fill(px,py,ii)
    end
    fplot(F,[xL(1),xL(2)]); xlabel('x'); ylabel('F(x)');
    disp(['area calculated by trapz.m for ',num2str(np),' points =',num2str(a2)]);
    title(['area calculated by trapz.m for ',num2str(np),' points =',num2str(a2)]);
    pause(dur);   % wait a bit
end

% ***************
% Method C
a1 = quad(F,xL(1),xL(2)); % use quad to integrate
msg = ['area calculated by quad.m = ' num2str(a1,10)]; disp(msg);