% exemplo de interpolação polinomial
% ARS {21/10/2005]

dx=0.04; % ruim para Ci*x^i
%dx=0.02; % ruim para extrapolação de Lagrange
%dx=0.01; % ruim para extrapolação de Newton e Neville
%dx=0.1;

x=[0.1:dx:0.9]'; % pontos nodais

phi=5;
y=sinh(phi*x)./(x*sinh(phi));

n=length(x);

xc=[0:0.01:1]'; % pontos para interpolação
m=length(xc);
yc(2:m)=sinh(phi*xc(2:m))./(xc(2:m)*sinh(phi));
yc(1)=phi/sinh(phi);

% interpolação polinomial por monômios
tic;
On=ones(n,1);
M=On;
for i=1:n-1
   M=[M x.^i];
end

C=inv(M)*y;  %coeficientes polinomiais com erros de arredondamento (sem pivotamento)
%C=M\y;  %coeficientes polinomiais sem erros de arredondamento (com pivotamento parcial)

yn1=C'*(On*xc').^([0:n-1]'*ones(1,m));

%for i=1:m
%   yn1(i)=C(n);
%   for j=n-1:-1:1
%      yn1(i)=yn1(i)*xc(i)+C(j);
%   end
%end

tp=toc

condM=cond(M)

% interpolação de Newton
% Tabela de Newton
tic;
A=zeros(n,n);
A(:,1) = y;
for i=2:n
   for j=1:n-i+1
      A(j,i)=(A(j+1,i-1)-A(j,i-1))/(x(i+j-1)-x(j));
   end
end

for i=1:m
   p=1;
   yn2(i)=A(1,1);
	for j=2:n
      p = p*(xc(i)-x(j-1));
      yn2(i) = yn2(i) + p * A(1,j);
	end
end
tn=toc

% interpolação de Lagrange
tic;
In=eye(n);
On=ones(1,n);

z=x*On;
dz=z'-z-diag(x);

for i=1:m
   w=z+xc(i)*In;
   yl1(i)=prod((xc(i)-w)./dz)*y;     %interpoladores de Lagrange
end
tl=toc

% interpolação de Neville
tic;
Q=zeros(n,n);
Q(:,1) = y;
for k=1:m
   for i=2:n
		for j=2:i
      	Q(i,j) = ((xc(k)-x(i-j+1))*Q(i,j-1)-(xc(k)-x(i))*Q(i-1,j-1))/(x(i)-x(i-j+1));
   	end
   end
   yn3(k) = Q(n,n);
end
tnl=toc

plot(xc,yc,'b',xc,yn1,'r',xc,yn2,'m',xc,yl1,'g',xc,yn3,'c',x,y,'o');
legend('exato','polinômio','Newton','Lagrange','Naville','raízes');
