%Comparing various numerical integration techniques, By Ryan Goh
%note: for a nice coding exercise, try to "vectorize" the code!!
 

%%Set function f, exact integral I, and interval x\in[a,b] to compare
f = @(xx) 1/xx; a = 1;b=4; I = log(4)
 %f = @(xx) exp(cos(xx)); a =0; b=1;  I = 2.341574841713053164922975186
%f = @(xx) cos(xx.^2); a = 0; b = 1; I = 0.904524237900272;

 
N = 2.^[2:1:12]';
ln = length(N);
%%%Collecting Errors
IM = zeros(ln,1);
IT = IM;
IS = IM;


EM = zeros(ln,1);
ET = zeros(ln,1);
ES = zeros(ln,1);

REM = EM;
RET = ET;
RES = ES;
 

for ii = 1:ln


%%Choose number of points
n = N(ii);  %For simpson's rule always choose n even!!!



%width
dx = (b-a)/n;

%%vector of partition points
x = [a:dx:b];

%%vector of midpoints
m = [a+dx/2:dx:b-dx/2];

S_mid=0;

%Calculate approximations

%%%%%%%
%%%%%%%%%midpoint
for jj=1:n
    S_mid = S_mid + f(m(jj));
end
I_mid = S_mid*dx


%%Error
AbsError_mid = abs(I_mid - I)

RelativeError_mid = abs(I_mid - I)/I

IM(ii) = I_mid;
EM(ii) = AbsError_mid;
REM(ii) = RelativeError_mid;

%%%%%%%%
%%%%%%%Trapezoid

S_trapz=f(x(1))/2 + f(x(n+1))/2;
for jj = 2:n
    S_trapz = S_trapz + f(x(jj));
end
I_trapz = S_trapz*dx
%%Error
AbsError_trapz = abs(I_trapz - I)

RelativeError_trapz = abs(I_trapz - I)/I

IT(ii) = I_trapz;
ET(ii) = AbsError_trapz;
RET(ii) = RelativeError_trapz;


%%%%%%
%%%%%%Simpson's Rule
S_Simp = f(x(1)) + f(x(n+1));
for jj = 2:2:n
    S_Simp = S_Simp + 4*f(x(jj));
end

for jj = 3:2:n-1
    S_Simp = S_Simp + 2*f(x(jj));
end

I_Simp = S_Simp*(dx/3)
AbsError_Simp= abs(I_Simp - I)

RelativeError_Simp = abs(I_Simp - I)/I


IS(ii) = I_Simp;
ES(ii) = AbsError_Simp;
RES(ii) = RelativeError_Simp;
end

figure(1)
plot(N,IM,'-o',N,IT,'-o',N,IS,'-o',N,I,'-o')
legend(['Midpoint'],['Trapezoid'],['Simpsons Rule'])

figure(2)
plot(N,EM,'-o',N,ET,'-o',N,ES,'-o')
title(['Absolute Error'])
xlabel('n')
ylabel('I - I_n')
legend(['Midpoint'],['Trapezoid'],['Simpsons Rule'])

figure(3)
plot(N,REM,'-o',N,RET,'-o',N,RES,'-o')
title(['Relative Error'])
xlabel('n')
ylabel('I - I_n')
legend(['Midpoint'],['Trapezoid'],['Simpsons Rule'])




LN = log(N);
figure(5)
plot(LN,log(EM),'-o',LN,log(ET),'-o',LN,log(ES),'-o')
title(['Absolute Error'])
xlabel('n')
ylabel('I - I_n')
legend(['Midpoint'],['Trapezoid'],['Simpsons Rule'])

figure(6)
plot(LN,log(REM),'-o',LN,log(RET),'-o',LN,log(RES),'-o')
title(['Relative Error'])
xlabel('n')
ylabel('I - I_n')
legend(['Midpoint'],['Trapezoid'],['Simpsons Rule'])




