
%interpolacja lagrange'a

%czy polyfit dziala

graphics_toolkit("gnuplot");
M=18;
x=1:5;
y=[-2 1 -5 3 -3]
p=polyfit(x,y,4);
s=linspace(1-0.1,5+0.1);
plot(x,y,"+",s,polyval(p,s));pause(1);

%zbieznosc dla sin i sin(4x) na wezlach rownoodleglych
a=-1;
b=4;
f1=@sin;
f2=@(x) sin(4*x);
k=1;
for n=6:4:M,
   x=linspace(a,b,n+1);
   y1=f1(x);
   p1=polyfit(x,y1,n);
   y2=f2(x);
   p2=polyfit(x,y2,n);
   s=linspace(a,b,10000);
   opis=[";wr n=", num2str(n), ";"]; 
   plot(x,y1,"+",s,p1s=polyval(p1,s),opis,s,f1s=f1(s),";f1;");pause(1);
   plot(x,y2,"+",s,p2s=polyval(p2,s),opis,s,f2s=f2(s),";f2;");pause(1);
   err(k++,:)=[n,norm(p1s-f1s,'inf'),norm(p2s-f2s,'inf')];
 end
 err
 
 %zbieznosc dla sin i sin(4x) na wezlach Czebyszewa
a=-1;
b=4;
f1=@sin;
f2=@(x) sin(4*x);
k=1;
for n=6:4:M,
   x=(a+b)/2 + 0.5*(b-a)*cos(pi*(0.5+(0:n))/(n+1));
   y1=f1(x);
   p1=polyfit(x,y1,n);
   y2=f2(x);
   p2=polyfit(x,y2,n);
   s=linspace(a,b,10000);
   opis=[";wc n=", num2str(n), ";"]; 
   plot(x,y1,"+",s,p1s=polyval(p1,s),opis,s,f1s=f1(s),";f1;");pause(1);
   plot(x,y2,"+",s,p2s=polyval(p2,s),opis,s,f2s=f2(s),";f2;");pause(1);
   erc(k++,:)=[n,norm(p1s-f1s,'inf'),norm(p2s-f2s,'inf')];
 end
 erc
 
 %wezly roownoodlegle na [-1,1.5] a blad interpolacji sin na [-1,4]
 
%zbieznosc dla sin i sin(4x) na wezlach rownoodleglych
a=-1;
b=4;
f=@sin;
k=1;
erh=[];
for n=6:4:M,
   x=linspace(a,0.5*(a+b),n+1);
   y=f(x);
   p=polyfit(x,y,n);
   s=linspace(a,b,10000);
   ps=polyval(p,s);
   fs=f(s);
   opis=[";wh n=", num2str(n), ";"]; 
   plot(x,y,"+",s,ps,opis,s,fs,";f1;");pause(1);
   erh(k++,:)=[n,norm(ps-fs,'inf')];
 end
 erh
 %nie widac zbieznosci ktora teoretycznie zachodzi....
 
 %przyklad Rungego
  %zbieznosc dla 1/(1+x^2) na [-5,5] na wezlach Czebyszewa
  %a rozbieznosc na rownoodleglych
a=-5;
b=5;
f1=@(x) 1.0./(1+x.*x);
k=1;
for n=2:4:M,
   x1=linspace(a,b,n+1);
   x2=(a+b)/2 + 0.5*(b-a)*cos(pi*(0.5+(0:n))/(n+1));
   y1=f1(x1);
   p1=polyfit(x1,y1,n);
   y2=f1(x2);
   p2=polyfit(x2,y2,n);
   s=linspace(a,b,10000);
   opis=["; n=", num2str(n), ";"]; 
   plot(x1,y1,"+;r;",s,p1s=polyval(p1,s),opis,s,f1s=f1(s),";f1;");pause(1);
   plot(x2,y2,"+;c;",s,p2s=polyval(p2,s),opis,s,f1s,";f2;");pause(1);
   erR(k++,:)=[n,norm(p1s-f1s,'inf'),norm(p2s-f1s,'inf')];
 end
 erR
 
 