%% Polynomial/Newton-Cotes Integration Examples
% Examples using data from textbook Example 4.2 and 5.1 (although with
% extra point)

format longg;

%Analytical Expressions
dx=0.1;
x=3.1:dx:3.9;
f=1./x;
df=-1./x.^2;

%Plot...
figure(1);
plot(x,f,x,df);
title('Analytical');
legend('f(x)','f''(x)');
xlabel('x'); ylabel('f(x) or f''(x)');
z
%Integral from 3.1 to 3.9

fp=log(abs(3.9/3.1));
err=0;
display(['Function Integral F(x)=ln(abs(x)) from x=3.1 to 3.9,', char(10), ' F(x)=', num2str(fp,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

%For midpoint rule, note that we first evaluate at the point values, not
%halfway between... so, we actually integrate from 3.05 to 3.95.
fm=dx*sum(f);
err=abs(fm-fp);
display(['Wrong Midpoint Rule from x=3.1 to 3.9,', char(10), ' F(x)=', num2str(fm,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

%What we actually did...
fp2=log(abs(3.95/3.05));
err=0;
display(['Function Integral F(x)=ln(abs(x)) from x=3.05 to 3.95,', char(10), ' F(x)=', num2str(fp2,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

%For midpoint rule, note that we could evaluate at the actual midpoints...
%We presumably don't have this data... But can evaluate for dx/2 offset limits.
fm=dx*sum(1./(x(1:end-1)+dx/2));
err=abs(fm-fp);
display(['Corrected Midpoint Rule from x=3.1 to 3.9,', char(10), ' F(x)=', num2str(fm,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

%For trapezoid rule, 
ft=sum((f(1:end-1)+f(2:end)).*dx./2);
err=abs(ft-fp);
display(['Trapezoid Rule from x=3.1 to 3.9,', char(10), ' F(x)=', num2str(ft,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

%For Simpson's rule, 
%  SO -- Think: How to code this as a loop?
ft=(1/3)*dx*(f(1)+4*f(2)+2*f(3)+4*f(4)+2*f(5)+4*f(6)+2*f(7)+4*f(8)+f(9));
err=abs(ft-fp);
display(['Simpson 1/3 Rule from x=3.1 to 3.9,', char(10), ' F(x)=', num2str(ft,'%1.16e'), ' error=', num2str(err,'%1.16e')]);
