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

clear; clc;

format longg;

x=[3.35 3.40 3.50 3.60]'
f=[0.298507 0.294118 0.285714 0.277778]'

%Direct / exact fit polynomial 
A=[ones(size(x)) x x.^2 x.^3]
b=A\f

xx=3.44
fn=1/xx;
err=0;
display(['Function f(x)=1/x at x=3.44,', char(10), ' f(x)=', num2str(fn,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fc=polyval(flipud(b),xx);
err=abs(fc-fn);
display(['Cubic Interpolation at x=3.44,', char(10), ' f(x)=', num2str(fc,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fq=polyval(vander(x(1:3))\f(1:3),xx);
err=abs(fq-fn);
display(['Quadratic Interpolation at x=3.44,', char(10), ' f(x)=', num2str(fq,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fi=polyval(vander(x(2:3))\f(2:3),xx);
err=abs(fi-fn);
display(['Linear Interpolation at x=3.44,', char(10), ' f(x)=', num2str(fi,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

fp=log(abs(3.6/3.3));  %MATLAB uses "log" for "ln", "log10" for "log"
err=0;
display(['Function Integral F(x)=ln(abs(x)) from x=3.3 to 3.6,', char(10), ' F(x)=', num2str(fp,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fcp=polyval(polyint(flipud(b)'),3.6)-polyval(polyint(flipud(b)'),3.3);
err=abs(fcp-fp);
display(['Cubic Interpolation and integration from x=3.3 to 3.6,', char(10), ' F(x)=', num2str(fcp,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fqp=polyval(polyint((vander(x(1:3))\f(1:3))'),3.6)-polyval(polyint((vander(x(1:3))\f(1:3))'),3.3);
err=abs(fqp-fp);
display(['Quadratic Interpolation and integration from x=3.3 to 3.6,', char(10), ' F(x)=', num2str(fqp,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

fip=polyval(polyint((vander(x(2:3))\f(2:3))'),3.6)-polyval(polyint((vander(x(2:3))\f(2:3))'),3.3);
err=abs(fip-fp);
display(['Linear Interpolation and integration from x=3.3 to 3.6,', char(10), ' F(x)=', num2str(fip,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

