%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) sqrt(1+xx.^4); a = 0; b = 6; I = 73.15271966572626;
%f = @(xx) cos(xx.^2); a = 0; b = 1; I = 0.9045242379002720814747;

 
%%%Collecting Errors
 
%%Choose number of points
n = 6;  %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

%%%%%%%%
%%%%%%%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



%%%%%%
%%%%%%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




