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

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=-1/xx^2;
err=0;
display(['Function Derivative f''(x)=-1/x^2 at x=3.44,', char(10), ' f''(x)=', num2str(fp,'%1.16e'), ' error=', num2str(err,'%1.16e')]);

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

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

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

%Lagrange Polynomial

a=x(1);
b=x(2);
c=x(3);

P1=((xx-b)/(a-b))*f(1)+((xx-a)/(b-a))*f(2);
err=abs(P1-fn);
display(['Lagrange 1st Order Interpolation (a,b) at x=3.44,', char(10), ' f(x)=', num2str(P1,'%1.16e'), ' error=', num2str(err,'%1.16e')]);
P1p=(1/(a-b))*f(1)+(1/(b-a))*f(2);
err=abs(P1p-fp);
display(['Lagrange 1st Order Interpolation (a,b) at x=3.44,', char(10), ' f''(x)=', num2str(P1p,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

P1=((xx-c)/(b-c))*f(2)+((xx-b)/(c-b))*f(3);
err=abs(P1-fn);
display(['Lagrange 1st Order Interpolation (b,c) at x=3.44,', char(10), ' f(x)=', num2str(P1,'%1.16e'), ' error=', num2str(err,'%1.16e')]);
P1p=(1/(b-c))*f(2)+(1/(c-b))*f(3);
err=abs(P1p-fp);
display(['Lagrange 1st Order Interpolation (b,c) at x=3.44,', char(10), ' f''(x)=', num2str(P1p,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

P2=((xx-b)*(xx-c)/((a-b)*(a-c)))*f(1)+((xx-a)*(xx-c)/((b-a)*(b-c)))*f(2)+((xx-a)*(xx-b)/((c-a)*(c-b)))*f(3);
err=abs(P2-fn);
display(['Lagrange 2nd Order Interpolation at x=3.44,', char(10), ' f(x)=', num2str(P2,'%1.16e'), ' error=', num2str(err,'%1.16e')]);
P2p=((2*xx-b-c)/((a-b)*(a-c)))*f(1)+((2*xx-a-c)/((b-a)*(b-c)))*f(2)+((2*xx-b-a)/((c-a)*(c-b)))*f(3);
err=abs(P2p-fp);
display(['Lagrange 2nd Order Interpolation at x=3.44,', char(10), ' f''(x)=', num2str(P2p,'%1.16e'), ' error=', num2str(err,'%1.16e'), char(10)]);

