

graphics_toolkit("gnuplot");

x=[-1,0,1,2,4]
y=[3,-1,1,2,10]
t=linspace(-1.1,4.1,300);

p=polyfit(x,y,length(x)-1);
plot(x,y,"*",t,polyval(p,t));pause(2);
y-polyval(p,x)

%splajny kubiczne not-a-knot

p1=spline(x,y);
y-ppval(p1,x)
plot(x,y,"*",t,polyval(p,t),";wiel int L ;",t,ppval(p1,t),";spline;");pause(2);

f=@sin;
p=polyfit(x,y=f(x),length(x)-1);
plot(x,y,"*",t,polyval(p,t),t,f(t),"b;sin;");pause(2);
y-polyval(p,x)

p1=spline(x,y);
y-ppval(p1,x)
plot(x,y,"*",t,ppval(p1,t),";spline;",t,f(t),"b;sin;");pause(2);


f=@(x) sin(4*x);
p=polyfit(x,y=f(x),length(x)-1);
plot(x,y,"*",t,polyval(p,t),t,f(t),"b;sin(4x);");pause(2);
y-polyval(p,x)

p1=spline(x,y);
y-ppval(p1,x)
plot(x,y,"*",t,ppval(p1,t),";spline;",t,f(t),"b;sin(4x);");pause(2);

%wielomian bazowy Lagrange'a l3 dla wezlow 0,...,10
%i splajn kubiczny s3zerujacy sie we wszystkich wezlach oprocz 3ciego 
%gdzie jest rowny 1 
x=0:10;
t=linspace(min(x),max(x),300);
y=zeros(size(x));
y(4)=1;
p=polyfit(x,y,length(x)-1);
p1=spline(x,y);
plot(x,y,"*",t,polyval(p,t),"r;l3 wezly rownoo;",t,ppval(p1,t),"b;spline s3;");
pause(2);

%wezly czebyszewa na [0,10]
a=0;b=10;n=10;
x=0.5*(a+b)-0.5*(b-a)*cos((0.5+(0:n))*pi/(n+1))
p=polyfit(x,y,length(x)-1);
p1=spline(x,y);
plot(x,y,"*",t,polyval(p,t),"r;l3 wezly czeb;",t,ppval(p1,t),"b;spline s3;");
pause(2);

%zbieznosc interpolacji Lagrange'a dla wezlow rownoodleglych Czebyszewa 



function er=testint(f,a,b,M=10,pplot=1e-2)
%function testing convergence of Lagrange interpolation
%at equidistant nodes and Chebyshev ones
%INPUT
%f- function handle to f - the interpolated function
%a,b - the interval endswitch
%M max degreee of polynomial (n+1 max no of nodes) default =10
%pplot - a threshold - when both discrete max errors are below pplot - there is no plot
%OUTPUT
%er - the matrix M x 3 - 1st column the degree of polynomial ; second col: the discrete max norm of the error 
%     between f and Lagrange polynomial on equid. nodes, 3rfd col: the discrete max norm of the error on the Chebyshev nodes
%EXAMPLE:
%  ers=testint(@sin,-1,4,4)
 x=linspace(a,b,1e4);
 fx=f(x);
 er=zeros(M,3);
 for n=1:M,
  er(n,1)=n; 
  wr=linspace(a,b,n+1);
  fr=f(wr);
  pr=polyfit(wr,fr,n);
  prx=polyval(pr,x);
  er(n,2)=norm(fx-prx,"inf");
  wc=0.5*(a+b)-0.5*(b-a)*cos((0.5+(0:n))*pi/(n+1));
  fc=f(wc);
  pc=polyfit(wc,fc,n);
  pcx=polyval(pc,x);
  er(n,3)=norm(fx-pcx,"inf");
  %[n,err=norm(er(n,2:end)),pplot]
  if(err=norm(er(n,2:end))>pplot)
   clf;
   hold on;
   plot(wr,fr,["*;n=",num2str(n),";"]);
   plot(x,fx,";f;",x,prx,";wr;");
   pause(1);
   clf;hold on;
   plot(wc,fc,["*;n=",num2str(n),";"]);
   plot(x,fx,";f;",x,pcx,";wc;");
   pause(2);
   hold off;
  else 
   disp("Norma ponizej progu");
  endif 
 endfor 
endfunction

ers=testint(@sin,-1,4,4);
ers3=testint(@(x) sin(3*x),-1,4);
erR=testint(@(x) 1./(1+x.*x),-5,5,18)




function [er,rat]=testspline(f,a,b,M=6,pplot=1e-3)
%function testing convergence of cubis spline interpolation on 
%at equidistant nodes (not-a -knot bnd conditions)
%INPUT
%f- function handle to f - the interpolated function
%a,b - the interval endswitch
%M max no of doubling the mesh
%pplot - a threshold - when both discrete max errors are below pplot - there is no plot
%OUTPUT
%er - the matrix M x 2 - 1st column the degree of polynomial ; 
%      second col: the discrete max norm of the error  between f and spline, 
%EXAMPLE:
%  ers=testspline(@sin,-1,4)
%rat - a vector with ratio of errors err(n,2)/err(n+1,2) 
 x=linspace(a,b,1e4);
 fx=f(x);
 er=zeros(M,2);
 N=8;
 for n=1:M,
  er(n,1)=n; 
  wr=linspace(a,b,N);
  fr=f(wr);
  pr=spline(wr,fr);
  prx=ppval(pr,x);
  er(n,2)=norm(fx-prx,"inf");
  %[n,err=norm(er(n,2:end)),pplot]
  if(er(n,2)>pplot)
   plot(wr,fr,["*;spline no of knots=",num2str(N),";"],x,fx,";f;",x,prx,";wr;");
   pause(2);
  else 
   %disp("Norma ponizej progu");
  endif 
  N=1+2*(N-1);
 endfor 
 rat=er(1:rows(er)-1,2)./er(2:rows(er),2); 
endfunction
[ers,rats]=testspline(@sin,-1,4)
printf("rzad zbieznosci dla sin(x):O(h^%g)\n",log2(rats(end)))
[ers3,rat3]=testspline(@(x) sin(3*x),-1,4)
printf("rzad zbieznosci dla sin(3x):O(h^%g)\n",log2(rat3(end)))
[erR,ratR]=testspline(@(x) 1./(1+x.*x),-5,5)
printf("rzad zbieznosci dla 1/(1+x^2):O(h^%g)\n",log2(ratR(end)))