

%interpolacja Lagrange'a

%test polyval narysuj wykres x^2-x na [-2,2] korzystajac z polyval
s=linspace(-2,2);
plot(s,polyval([1,-1,0],s),";wykres x^2-x;");pause(2);


%znajdz wiel interp Lagrange'a dla tabelki
%x(k): 0  1 2  4
%y(k)  -1 7 4 -2 
%narysuj wykres zaznaczjac pkty interpolacji (x(k),y(k) gwiazdka
x=[0 1 2 4];
y=[-1 7 4 -2];
a=polyfit(x,y,length(x)-1)
s=linspace(min(x)-1,max(x)+1);
plot(x,y,"*",s,polyval(a,s),";int L na 4p;");pause(2);
%to samo ale zamien wezel 4 z 2.1
x(4)=2.1;

a=polyfit(x,y,length(x)-1)
s=linspace(min(x)-1,max(x)+1);
plot(x,y,"*",s,polyval(a,s));pause(2);

%dla wezlow -1,0,1 lub -1,0,1,10
%narysuj wykres wiel interp Lagrange'a w dla f(x)=sin(x)
%policz przyblizonya norme max f-w na [min(x),max(x)]
a=-1;b=10;
x=[-1,0,1,10];
y=sin(x);
p=polyfit(x,y,length(x)-1);
s=linspace(a,b,1000);
ys=sin(s);
ls=polyval(p,s);
er=norm(ys-ls,"inf")
plot(x,y,"+",s,ys,s,ls,";ls;");pause(2);

%narysuj dla wezlow -1,2,3,4 f bazowa Lagrange'a
% dla wezla 3
x=[-1 2 3 4]
y=zeros(size(x))
y(3)=1;
p=polyfit(x,y,length(x)-1);
s=linspace(-1,4);
plot(x,y,"b*",s,polyval(p,s),";f baza lagrange'a l_3;");pause(2);

%przetestowac czy ciag wiel Lagrange'a na wezlach ronoodleglych 
%zbiega na [0,2pi] dla sin(x) i sin(4x)

function err=testbladiL(f,a,b,typ=0,M=20,disp=-1)
%testujemy blad interp Lagrange'a  dal kolejnej ilosci wezlow od 2 do M+1
%INPUT
%f - uchwyt do funkcji
%a,b - konce odcinka 
%typ - typ wezlow 0 rownoodlegle 1 Czebyszewa
%M - max stopien wielomianu interpolacyjnego
%disp - max stopien wielomiany 
%        dla ktorego rysowane sa wykresy f i wiel interp a wykres bledu 
n=2;
err=zeros(M-1,3);
s=linspace(a,b,1e5);
fs=f(s);
for k=2:M,
  if(typ==0)
    x=linspace(a,b,n+1);
  else 
    x= 0.5*(a+b) + 0.5*(b-a)*cos(pi*(0.5+(0:n))/(n+1));  
  end  
 y=f(x);
 p=polyfit(x,y,n);
 ps=polyval(p,s);
 err(k-1,1)=n;
 rei=err(k-1,2)=norm(polyval(p,x)-y,"inf");
 er=err(k-1,3)=norm(fs-ps,"inf");
 ers=fs-ps;
 if(n<disp)
   plot(s(1:100:end),(fs-ps)(1:100:end),[";error n=",num2str(n)," typ=",num2str(typ),";"],x,0*x,"g*");pause(1);
   plot(s(1:100:end),fs(1:100:end),";f;",s,ps,[";n=",num2str(n)," typ=",num2str(typ),";"],x,f(x),"g*");pause(1);
 end  
 %err(k-1,:)
 n++;
end 
end


f=@sin;
a=0;b=2*pi;M=4;typ=0;disp=5;
err=testbladiL(f,a,b,typ=0,M=20,5)
f=@sin;
a=0;b=2*pi;
erc=testbladiL(f,a,b,typ=1,M=20,5)

%przyklad Rungego

f=@(x) 1.0./(1+x.*x);
a=-5;b=5;
err=testbladiL(f,a,b,typ=0,M=25,5)
erc=testbladiL(f,a,b,typ=1,M=25,5)


%interpolacja hermite n wezlow 2 krotnych
function [p,A]=iHw2kr(x,y,dy,typ=0)
%f znajduje wspolczynniki wiel interp Hermite'a w w bazie (x^(N),...,x,1)
%dla wezlow w x takich że w(x(k)=y(k) w'(x(k)=dy(k) k=1:length(x) type=0
%lub w(x(k)=y(k)  k=1:length(x) type=0 w'(x(k)=dy(k) k=2:length(x)-1
%czyli wewnetrzne (o ile monotoniczne...) wezly sa 2krote, koncowe 1 krotne 
%INPUT%
%x- wektor (kolumna) z welzami
%y - wektor z wartosciami
%dy - wektor z wartosciami pochodnej w wezlach
%OUPUT
%p - wspolczynniki wielomianu w bazie (x^(2n-1),...,x,1)
%A - macierz pozwalajaca znalezc wspolczynniki

n=length(x);
A=fliplr(vander([x;x]));
A2=zeros(n,2*n);
A2(:,2:2*n)=A(1:n,1:(2*n-1))*diag(1:(2*n-1));
A((n+1):2*n,:)=A2;
A=fliplr(A);

if(typ==0)
 p=A\[y;dy]
else
  N=length(A);
  A=A(2:N-1,2:N-1);
  %p=A\[y;dy(2:n-1)];
end 
end



a=0;b=2;
x=[a,0.5*(a+b),b]'
f=@(x) -sin(x);
df=@(x) -cos(x);
y=f(x);
dy=df(x);
p=iHw2kr(x,y,dy);
xx=linspace(a-1,b+1,300);

plot(xx,f(xx),";f;",xx,polyval(p,xx),";wiel int Herm;",x,y,"g+");
pause(1);
