%%Numerical Integration with integral

% Integration is one of the two fundamental operations in calculus
% If you aren't familiar with calculus, we can think of integration as
% finding the "area under a curve"
% Calculates the area from a given starting point to a given end point

%Built in function in MATLAB - integral(fun, xmin, xmax)
%fun is a function handle to the function that we want to integrate (find
%area under)
%We can use anonymous functions for fun, allows us to avoid creating new
%files for each function
%Note: fun must use elementwise operations

%xmin is the start of the region we are integrating over
%xmax is the end of the region we are integrating over
%Note that xmin can take the value -Inf (negative infinity) and xmax can
%take the value Inf (positive infinity) which allows you to calculate an
%improper integral

%Two curves we will use for integration
xx = 0 : 0.1 : 5;
yy1 = (xx - 3).^2 + 5;
yy2 = xx.*sin(xx) + 10;
plot(xx, yy1, 'b',  xx, yy2, 'r')
axis([0 5 0 15])

%Example 1: find area under y1 between 1 and 4
y1 = @(x)(x-3).^2 + 5;
area1 = integral(y1,1,4);
disp(area1)

%Example 2: find area between y1 and y2
%Need to find endpoints of the region
%Can do this by finding where y1 = y2 -> y1 - y2 = 0
%So we can use fzero!
y2 = @(x)(x.*sin(x)+10);
difffunc = @(x)y1(x)-y2(x);
root1 = fzero(difffunc,0.75);
root2 = fzero(difffunc,4.25);
%disp(root1)
%disp(root2)
%Find area between the curves by taking the difference -> area under the
%top curve - area under the bottom curve
areabetween = integral(y2,root1,root2) - integral(y1,root1,root2);
disp(areabetween)


%Example 3: find b such that the area under y2 between 0 and b is equal to
%30
%Step 1: Define a function that calculates the area under the curve from 0
%to b
area0tob = @(b)integral(y2,0,b);
%Step 2: To find (area 0 to b) = 30, can use fzero to calculate (area 0 to b) - 30 =
%0
%Step 3: Using fzero, need a starting point. Plot area0tob-30.
x = 0:.1:5;
y = zeros(length(x));
for i=1:length(x)
  y(i) = area0tob(x(i))-30;
end
plot(x,y)
%2.5 seems like a good starting point
result = fzero(@(x)area0tob(x)-30,2.5);
disp(result)