%% Polynomial Approximations
% Demonstration of difference based polynomials

%% Neville's Algorithm Demo

% x=[3.20 3.30 3.35 3.40 3.50 3.60 3.65 3.70]';
% f=[0.312500 0.303030 0.298507 0.294118 0.285714 0.277778 0.273973 0.270270]';
x=[3.40 3.50 3.35 3.60]';
f=[0.294118 0.285714 0.298507 0.277778]';
N=size(x,1);
X=3.44;

g=zeros(N);
g(:,1)=f;
format longg;

for n=1:N
    for i=1:(N-n)
        g(i,n+1)=(((X-x(i))*(g(i+1,n))-((X-x(i+n))*(g(i,n))))./(x(i+n)-x(i)));
    end;
end;

g

%% Forward Difference Demo
x=[3.20 3.30 3.40 3.50 3.60 3.70]';
f=[0.312500 0.303030 0.294118 0.285714 0.277778 0.270270]';

N=size(f,1);

df=zeros(N);
df(:,1)=f;

for n=1:N
    for i=1:(N-n)
        df(i,n+1)=df(i+1,n)-df(i,n);
    end;
end;

df

s=(X-x(1))./(x(2)-x(1));
%Note: Need to use definition in terms of gamma functions n!=gamma(n+1)
C1=gamma(s+1)./(gamma(1+1)*gamma(s-1+1));
C2=gamma(s+1)./(gamma(2+1)*gamma(s-2+1));
C3=gamma(s+1)./(gamma(3+1)*gamma(s-3+1));

p=df(1,1)+C1*df(1,2)+C2*df(1,3)+C3*df(1,4)

%% Spline Demo

PS=spline(x,f)
figure(1);
xx=3:0.01:4;
plot(xx,ppval(PS,xx)); hold on; scatter(x,f);
