% This script calculates the area under a user defined curve, with user
% defined start and endpoints and stepsize
% f is the polynomial function that the user inputs
% h is the step size for the rectangles
% x1 and x2 are the limits on the graph that the user wishes to integrate
% over.
%
% Essentially, this script finds the antiderivative of the polynomial function using
% polyfit after finding the integral with predefined limits, it plots the
% original function and the integrated function, and it calculates the area
% of the integral.

% All of the functions and loops I used were smashed into this one Mega-script

clear all
clf % Clear the figure if any are open, just in case.
syms x

% % user-inputs if wanted
% f = input('Enter f(x)(polynomial only): ');
% h = input('Enter stepsize (h) of the function: ');
% x1 = input('Left limit of graph: ');
% x2 = input('Right limit of graph: ');

%premade inputs with polynomial
f = x^2 + 5*x + 10 %this example shows different degrees of antidifferentiation
h = .005; % a larger step size will increase the time the script takes to run, but provide much more accurate results.
x1 = -10; %from -10 to 10, nothing too fancy like a -1 two +3
x2 = 10;
%Plugged into the calculator, an answer of 8.6666e+02 was found compared to
%this script's best answer of 8.6746e+02 - rather close for an approximation

% f0 = 10; was going to try to use an initial condition to solve for C, but
%no luck

% Telling matlab that functions are functions.
f = @(x) eval(f);

hold on

rect = (x2 - x1)/h + 1;
a = x1; % Starting point in the integration.
x = x1:h:x2;


% ----------------------------------------------------------
% Function/loop #1
% Using the Midpoint Rule to find the antiderivative
% This loop should solve the integral using a sum of small rectangles,
% aka the midpoint rule
% M = (b-a)*f((a+b)/2);

% The loop that calculates the area.
for ii = 1:1:rect
    b = a + h;
    M(ii) = ((b-a)/2)*(f(a)+f(b));
    If(ii) = sum(M);
    a = b;
end

% % Kinda gave up on this, unfortunately; ran out of time and patience
% %trying to solve for C with an initial condition
% If(1000)  %[(x2 - x1)/h] / 2
% c = If(1000) - f0
% [minval, pos]=min()

plot(x,f(x),x,If) % Plots the graph
xlabel('X axis')
ylabel('Y axis')
title('Comparison of function and integrated function')
legend('Original function','Antiderivative function')
    S = sum(M);
    
n = x(1:1:rect); % Establishing the vector lengths.
m = If(1:1:rect);

clear x % Establishing x as a variable again.
syms x
syms c

% Using the polyfit function, we can have Matlab find the function of the
% antiderivative using points found from the integration loop.

% A much more accurate antiderivative is found by using a smaller step
% size, at the cost of time that the script takes to run.
p = @(x) eval(p);
p = polyfit(n,m,6); % The six is for 6 degrees of polynomials, larger polynomials will need a larger buffer.
digits(3)
for ii = 1:length(p) % removes the smaller values of x in the antidifferentiation
   if p(ii) < 1e-10
       p(ii) = 0;
   end
end
display('The antiderivative of your function is:')
vpa(p(1),2)*x.^6 + vpa(p(2),2)*x.^5 + vpa(p(3),2)*x.^4 + vpa(p(4),2)*x.^3 + vpa(p(5),2)*x.^2 + vpa(p(6),2)*x + c
% vpa allows for some degree of precision when finding the antiderivative,
% without it, antiderivatives get ugly fast.
% This displays the antiderivative for the user.


% Now onto the wonderful world of numerical integration
    
% ----------------------------------------------------------
% Function Loop #2
% Using Midpoint Rule to find area
% This loop should solve the integral using a sum of small rectangles
% aka the midpoint rule
% M = (b-a)*f((a+b)/2);

% Have to clear everything from the previous loop to start from fresh

clear ii
clear b
clear M(ii)
clear If(ii)
clear S

rect = (x2 - x1)/h;
a = x1;
x = x1:h:x2;


% The loop that calculates the area.
for ii = 1:1:rect
    b = a + h;
    M(ii) = ((b-a)*f((a+b)/2));
    If(ii) = sum(M);
    a = b;
end

     S = sum(M);
    sprintf('Your area under your function and between your endpoints using the midpoint rule is %d\n',S)


%-----------------------------------------------------------
% Function Loop #3
% Trapezoidal Rule
% Trapezoidal rule makes a bunch of trapezoids along the curve, and then
% uses that to approximate the area.

% Have to clear everything from the previous loop to start from fresh.
clear ii
clear b
clear M(ii)
clear If(ii)
clear S

rect = (x2 - x1)/h; % Establishing step size and initial points again.
a = x1;
x = x1:h:x2;

% The loop that calculates the area with the trapezoidal method.
for ii = 1:1:rect
    b = a + h;
    M(ii) = ((b-a)/2)*(f(a)+f(b));
    S(ii) = sum(M);
    a = b;
end
S = sum(M);
    sprintf('Your area under your function and between your endpoints using the trapezoidal rule is %s\n',S)

   
% ----------------------------------------------------------
% Function Loop #4
% Simpson's Rule
% Simpson's rule makes a really close function to mimic the user inputted
% function, and then finds the area under the new function made to
% approximate the inputted function. Somewhat resembles a Taylor series.

clear ii % Refreshing the loop.
clear b
clear M(ii)
clear If(ii)
clear S

rect = (x2 - x1)/h;
a = x1;
x = x1:h:x2;
% The loop that calculates the area.

for ii = 1:1:rect
    b = a + h;
    M(ii) = (b-a)/6*(f(a) + 4*f((a+b)/2) + f(b));
    S(ii) = sum(M);
    a = b;
end
S = sum(M);
    sprintf('Your area under your function and between your endpoints using Simpsons Rule is %s\n',S)
    
% Nothing else beyond here, though.