% NUMERICAL INTEGRATION EXAMPLES WITH MATLAB
%   (Code by Ben Walter, METU-NCC, benjamin@metu.edu.tr)
%
% CONTENTS:
%
% Integrating f(x)
% -- via the right-endpoint rule
% -- via the left-endpoint rule
% -- via the midpoint rule
% -- via the trapezoidal rule
% -- via simpson's rule
% -- comparison of methods (Optional)
%
% Line Integrals of Vector Fields
%
% Line Integrals of Functions (Optional)
%
% Integrating f(x,y)
% -- on a rectangular region
% -- on a general region
%
% Example: Compute Center of Mass
%
% Example: Compute Area via Green's Theorem
%
%% INTEGRATING F(X) VIA RIGHT-ENDPOINTS
%
clc

disp('We will approximate the integral of f(x) = sin^2(x) from 1 to 3')
disp(' using 100 rectangles')

f = inline('(sin(x)).^2', 'x') % = function to integrate
a = 1                          % = starting point of integration
b = 3                          % = ending point of integration

n = 100                        % = number of rectangles used in approximation

%%
dx = (b-a)/n                   %  = width of each rectangle
right_endpoints = (a+dx):dx:b; %  = vector of right endpoints of rectangles
                               %     (format is  "<start>:<stepsize>:<end>")
                              
   % the command "sum(<vector>)" adds up all elements of a vector
right_approx_integral = sum( f(right_endpoints) )*dx   % int f(x) dx




%% INTEGRATING F(X) VIA LEFT-ENDPOINTS
%

dx = (b-a)/n;                  % = width of each rectangle
left_endpoints = a:dx:(b-dx);  % = vector of left endpoints of rectangles

left_approx_integral = sum( f(left_endpoints) )*dx     % int f(x) dx




%% INTEGRATING F(X) VIA MIDPOINTS

dx = (b-a)/n;                      % = width of each rectangle
midpoints = (a+dx/2):dx:(b-dx/2);  % = vector of midpoints of rectangles

midpoint_approx_integral = sum( f(midpoints) )*dx      % int f(x) dx




%% INTEGRATING F(X) VIA TRAPEZOIDS

dx = (b-a)/n;        % = width of each rectangle
endpoints = a:dx:b;  % = vector of endpoints of rectangles

                        % use ... to continue long formulas on the next line
trapezoid_approx_integral = ...
    sum( f(endpoints(2:end-1)) )*dx + ...
    f(endpoints(1))*dx/2 + ...   % instead of skipping one of the endpoints
    f(endpoints(end))*dx/2       % like in the right or left endpoint rule,
                                 % we'll keep them both, but only have
                                 % them each count half as much

% To be most careful, numerically, we should have added everything up before
% multiplying by dx -- that way we would reduce roundoff error.  To add
% everything up, we can use the code:
%    sum([ f(endpoints(1))/2, f(endpoints(2:end-1)), f(endpoints(end))/2 ])


%% INTEGRATING F(X) VIA SIMPSON

n = ceil(n/2)*2;     % we need n to be even for simpson's rule
dx = (b-a)/n;        % = width of each rectangle
endpoints = a:dx:b;  % = vector of endpoints of rectangles

simpson_approx = ...
    sum( f(endpoints(2:2:end-1)) )*4/3 + ...  %  even terms count 4/3
    sum( f(endpoints(3:2:end-2)) )*2/3 + ...  %  odd terms count 2/3
    f(endpoints(1))/3 + ...                   %  ends count 1/3
    f(endpoints(end))/3;                      %  ends count 1/3
simpson_approx = simpson_approx*dx



%% COMPARISON OF METHODS
clc
disp('We will compare the approximations given by different methods.')
disp('  (You should try this with different functions and limits of integration.)')

f = inline('(sin(x)).^2', 'x')
b = 3
a = 1

%%
disp('             Table of approximations given by different methods:')
disp('  n  | left-endpt | right-endpt | midpoint | trapezoid | simpson  |  random')
disp(' ----|------------|-------------|----------|-----------|----------|---------')

for (n=[2:4:40 60:60:800])  % estimate using 2, 6, 10, 14,... rectangles

    dx=(b-a)/n;                      % = width of rectangles
    endpoints=a:dx:b;                % = endpoints of rectangles
    midpoints=a+dx/2:dx:b-dx/2;% = midpoints of rectangles
    
    left_endpoint_approx = sum( f(endpoints(1:end-1)) )*dx;
    right_endpoint_approx = sum( f(endpoints(2:end)) )*dx;
    midpoint_approx = sum( f(midpoints) )*dx;
    trapezoid_approx = ...
        sum([ f(endpoints(1))/2, f(endpoints(2:end-1)), f(endpoints(end))/2 ] )*dx;
	simpson_approx = ...
        sum( f(endpoints(2:2:end-1)) )*4/3 + ...  %  even terms count 4/3
        sum( f(endpoints(3:2:end-2)) )*2/3 + ...  %  odd terms count 2/3
        f(endpoints(1))/3 + ...                   %  ends count 1/3
        f(endpoints(end))/3;                      %  ends count 1/3
    simpson_approx = simpson_approx*dx;
    
     % In theory we can also approximate the integral by just picking a
     % bunch of random points.....  let's try that, too!
    randompoints = a + (b-a).*rand(n,1);  % n random points between a and b
    random_approx = sum( f(randompoints) )*dx;
    
    fprintf('%3d  |  %f  |   %f  | %f |  %f | %f | %f\n', n, left_endpoint_approx, right_endpoint_approx, midpoint_approx, trapezoid_approx, simpson_approx, random_approx)
end

 % we can use the matlab command "quadl" to get matlab to compute an
 % approximation using its hard-coded black magic....
matlab_approx = quadl(f,1,3);
fprintf('\n Actual value: %f \n', matlab_approx)




%% LINE INTEGRALS OF VECTOR FIELDS

clc

disp('F = [-y x]  calculate line integral around circle r(t) = [cos(t) sin(t)]')

P = inline('-y', 'x', 'y')     % F = [P Q]  The vector field to be
Q = inline('x', 'x', 'y')      %              integrated

x = inline('cos(t)', 't')      % Parameterization of the curve
y = inline('sin(t)', 't')      %  we are integrating along

a = 0                          % Bounds for the parameter t
b = 2*pi

n = 5                          % Number of rectangles in t-coordinate

disp('First we''ll show everything happening:')

%%
dt = (b-a)/n                             % width in t-coordinate

endpoints = a:dt:b                       % endpoints in t-coordinate
midpoints = a+dt/2:dt:b-dt/2             % We will use the midpoint rule.

dx = x(endpoints(2:end)) - x(endpoints(1:end-1)) % widths in x-coordinate
dy = y(endpoints(2:end)) - y(endpoints(1:end-1)) % widths in y-coordinate

line_integral = sum([ P(x(midpoints), y(midpoints)).*dx ... % int P dx + Q dy
                      Q(x(midpoints), y(midpoints)).*dy ])


%%
disp('Now let''s compute seriously:')
n = 1000

%%
dt = (b-a)/n;                            % width in t-coordinate

endpoints = a:dt:b;                      % endpoints in t-coordinate
midpoints = a+dt/2:dt:b-dt/2;            % We will use the midpoint rule.

dx = x(endpoints(2:end)) - x(endpoints(1:end-1)); % widths in x-coordinate
dy = y(endpoints(2:end)) - y(endpoints(1:end-1)); % widths in y-coordinate

line_integral = sum([ P(x(midpoints), y(midpoints)).*dx ... % int P dx + Q dy
                      Q(x(midpoints), y(midpoints)).*dy ])


% NOTE: This gave circulation.  We could calculate flux by solving 
%   int Q dx - P dy
% instead.
                  
                  
%% LINE INTEGRALS OF FUNCTIONS (OPTIONAL)

clc

disp('f(x,y) = cos(x^2-y) calculate line integral around circle r(t) = [cos(t) sin(t)]')

f = inline('cos(x.^2-y)', 'x', 'y') % function to be integrated

x = inline('cos(t)', 't')           % Parameterization of the curve
y = inline('sin(t)', 't')           %  we are integrating along

a = 0                               % Bounds for the parameter t
b = 2*pi

n = 5                               % Number of rectangles in t-coordinate

disp('First we''ll show everything happening:')

%%
dt = (b-a)/n                             % width in t-coordinate

endpoints = a:dt:b                       % endpoints in t-coordinate
midpoints = a+dt/2:dt:b-dt/2             % We will use the midpoint rule.

dx = x(endpoints(2:end)) - x(endpoints(1:end-1)) % widths in x-coordinate
dy = y(endpoints(2:end)) - y(endpoints(1:end-1)) % widths in y-coordinate
ds = sqrt( dx.^2 + dy.^2 )                       % widths in arc length coordinate

line_integral = sum( f(x(midpoints), y(midpoints)).*ds )   % int f(x,y) ds



%%
disp('Now let''s compute seriously:')
n = 1000

%%
dt = (b-a)/n;                             % width in t-coordinate

endpoints = a:dt:b;                       % endpoints in t-coordinate
midpoints = a+dt/2:dt:b-dt/2;             % We will use the midpoint rule.

dx = x(endpoints(2:end)) - x(endpoints(1:end-1)); % widths in x-coordinate
dy = y(endpoints(2:end)) - y(endpoints(1:end-1)); % widths in y-coordinate
ds = sqrt( dx.^2 + dy.^2 );                       % widths in arc length coordinate

line_integral = sum( f(x(midpoints), y(midpoints)).*ds )   % int f(x,y) ds



                  

%% INTEGRATING F(X,Y) ON A RECTANGLE
%
clc

disp(' Estimate the integral of cos(x^2-y) over the rectangle [a,b]x[c,d]')

f = inline( 'cos(x.^2 - y)', 'x', 'y' )
a = 1
b = 4
c = -1
d = 3

n = 5   % number of divisions in x-direction
m = 6   % number of divisions in y-direction

%%
fprintf('-------\n First let''s show everything that is happening:')
dx = (b-a)/n    % = width in the x-direction
dy = (d-c)/m    % = width in the y-direction

x_midpoints = a+dx/2:dx:b-dx/2  % = x-coordinates of midpoints
y_midpoints = c+dy/2:dy:d-dy/2  % = y-coordinates of midpoints

 % the matlab command "meshgrid" converts vectors of x and y coordinates
 % into a matrix of all points with these coordinates
[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints)

 % since midpoints_x and midpoints_y are matrices, we must use the command
 % sum(sum(<..>)) to add up all elements -- the first sum adds up the
 % columns, and the next sum adds up the results
midpoint_approx = sum(sum(f(midpoints_x, midpoints_y)))*dx*dy % int int f dx dy

%%
fprintf('-------\n Now let''s get serious about approximating the integral:')
n=500
m=600

%%
dx = (b-a)/n;
dy = (d-c)/m;

x_midpoints = a+dx/2:dx:b-dx/2;
y_midpoints = c+dy/2:dy:d-dy/2;

[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints);

midpoint_approx = sum(sum(f(midpoints_x, midpoints_y)))*dx*dy


disp('Try this for different functions, limits, and numbers of divisions!')

%  NOTE:  We could also have approximated this double integral using the
%  trapezoid rule.  For double integrals on rectangular regions, the
%  trapezoid rule says to add up f(endpoints_x, endpoints_y) weighting the
%  corners by 1/4, the edges by 1/2 and the middle by 1.



%% INTEGRATING F(X,Y) ON A REGION
%
clc

disp(' Estimate the integral of xy over the upper-right quadrant of a circle')

f = inline( 'x.*y', 'x', 'y')       % function to integrate
a = 0              % We will use a rectangle surrounding our region of 
b = 1              % integration in order to build our mesh
c = 0              % The rectangle is: [a,b]x[c,d]
d = 1

n = 5
m = 3

g = inline( 'x.^2 + y.^2', 'x', 'y')  % boundary of the region

%%
fprintf('-------\n First let''s show everything that is happening:')
dx = (b-a)/n
dy = (d-c)/m

x_midpoints = a+dx/2:dx:b-dx/2
y_midpoints = c+dy/2:dy:d-dy/2

[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints)

interior_midpoints_x = midpoints_x(g(midpoints_x, midpoints_y) <= 1)
interior_midpoints_y = midpoints_y(g(midpoints_x, midpoints_y) <= 1)

midpoint_approx = sum(f(interior_midpoints_x, interior_midpoints_y))*dx*dy

%%
fprintf('-------\n Now let''s get serious about approximating the integral:')
n=500
m=600

%%
dx = (b-a)/n;
dy = (d-c)/m;

x_midpoints = a+dx/2:dx:b-dx/2;
y_midpoints = c+dy/2:dy:d-dy/2;

[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints);

interior_midpoints_x = midpoints_x(g(midpoints_x, midpoints_y) <= 1);
interior_midpoints_y = midpoints_y(g(midpoints_x, midpoints_y) <= 1);

midpoint_approx = sum(f(interior_midpoints_x, interior_midpoints_y))*dx*dy

disp('Try this for different functions, limits, and numbers of divisions!')



%% EXAMPLE:  COMPUTE THE CENTER OF MASS OF THE TOP-RIGHT QUARTER OF A CIRCLE

clc

disp(' Compute center of mass upper-right quadrant of a circle:')

moment_x = inline( 'x', 'x', 'y');    % = vector of x-values
moment_y = inline( 'y', 'x', 'y');    % = vector of y-values
density = inline( 'x-x+1', 'x', 'y'); % = vector with 1 for every value of x
g = inline( 'x.^2 + y.^2', 'x', 'y'); % boundary of the region

a=0; b=1;                        c=0; d=1;
n=1000;                          m=1000;
dx=(b-a)/n;                      dy=(d-c)/m;
x_midpoints = a+dx/2:dx:b-dx/2;  y_midpoints = c+dy/2:dy:d-dy/2;

[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints);

interior_midpoints_x = midpoints_x(g(midpoints_x, midpoints_y) <= 1);
interior_midpoints_y = midpoints_y(g(midpoints_x, midpoints_y) <= 1);

XMoment_approx = sum(moment_x(interior_midpoints_x, interior_midpoints_y))*dx*dy
YMoment_approx = sum(moment_y(interior_midpoints_x, interior_midpoints_y))*dx*dy
Mass_approx = sum(density(interior_midpoints_x, interior_midpoints_y))*dx*dy

fprintf(' Center of mass: (%f, %f)\n', XMoment_approx/Mass_approx, YMoment_approx/Mass_approx)




%% EXAMPLE:  USE GREEN'S THEOREM TO COMPUTE AREA
%
clear all;
clc

r = @(t) (1 + 0.1*cos(7*t))       % a curve written in polar form

disp(' The area inside the polar curve above is difficult to get using double integrals')
disp(' Instead we''ll use Green''s theorem with vector field F=[0 x] :')
disp('   int_C x dy = double int_R d/dx(x) dx dy')
disp('              = double int_R  1      dx dy')
disp('              = Area(R)')

a = 0; b = 2*pi;                  % limits of t
n = 1000;                         % number of segments

x = @(t) (r(t).*cos(t))           % parameterize curve in rectangluar coords
y = @(t) (r(t).*sin(t))

P = inline('0', 'x', 'y')         % F = [P Q]  the vector field 
Q = inline('x', 'x', 'y')         %             to be integrated

%%

dt = (b-a)/n;                            % width in t-coordinate

endpoints = a:dt:b;                      % endpoints in t-coordinate
midpoints = a+dt/2:dt:b-dt/2;            % We will use the midpoint rule.

dx = x(endpoints(2:end)) - x(endpoints(1:end-1)); % widths in x-coordinate
dy = y(endpoints(2:end)) - y(endpoints(1:end-1)); % widths in y-coordinate

line_integral = sum([ P(x(midpoints), y(midpoints)).*dx ... % int P dx + Q dy
                      Q(x(midpoints), y(midpoints)).*dy ])
                  
         % Note:  Since P is zero, we could have skipped the dx stuff.

%%
%  Now let's try to compute the area using a double integral

one = inline( 'x-x+1', 'x', 'y'); % = vector with 1 for every value of x

a=-1.2; b=1.2;                   c=-1.2; d=1.2;
n=1000;                          m=1000;     % We'll need lots of memory!
dx=(b-a)/n;                      dy=(d-c)/m;
x_midpoints = a+dx/2:dx:b-dx/2;  y_midpoints = c+dy/2:dy:d-dy/2;

[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints);

theta = @(x,y) (atan(y./x)+pi*(x<0));  % get angle for (x,y) coordinate

interior_midpoints_x = midpoints_x(sqrt(midpoints_x.^2+midpoints_y.^2) <= r(theta(midpoints_y,midpoints_x)));
interior_midpoints_y = midpoints_y(sqrt(midpoints_x.^2+midpoints_y.^2) <= r(theta(midpoints_y,midpoints_x)));

double_integral = sum(one(interior_midpoints_x, interior_midpoints_y))*dx*dy

disp('It took a bit more time for the computer to get it this way...')


%% SURFACE INTEGRALS OF FUNCTIONS (OPTIONAL)

f = inline('x*y^2 - z', 'x', 'y', 'z')

x = inline('u', 'u', 'v')               %  parameterization of surface
y = inline('u.^2', 'u', 'v')            %    r(u,v) = [u  u^2  v]
z = inline('v', 'u', 'v')               %  Kreyszig p 451

a = 0   %  lower limit of u
b = 2   %  upper limit of u
c = 0   %  lower limit of v
d = 3   %  upper limit of v

disp('Let''s show everything this first time...')

n = 5   %  divisions for u
m = 3   %  divisions for v

%%

du = (b-a)/n
dv = (d-c)/m

% .... (to be added)


%% EXAMPLE: FLUX ACROSS A SURFACE (OPTIONAL)


% .... (to be added)


%% EXAMPLE: DIVERGENCE THEOREM (OPTIONAL)


% .... (to be added)




