

#a stworz splajn dla wezlow -1,0,1,3 z wartosciami w tych wezlach 
#1,7,-3,2 (czy innymi) typu not-a-knot
#oraz hermitowski przyjmujac zerowe wartosci pochdnych w koncach
#narysuj wykresy nw 1 okienku

x=[-1,0,1,3];
y=[1,7,-3,2];
pp1=spline(x,y);
pp2=spline(x,[0,y,0]);
s=linspace(-1,3);
plot(s,ppval(pp1,s),s,ppval(pp2,s),x,y,"+")
pause(1);
#wykres blisko konca lewego
plot(s=s(1:10),ppval(pp1,s),s,ppval(pp2,s))
pause(1);

#zd2 splajn na wezlach -5,-4,..,0..,5 taki ze
#(a) s(1)=1 s(j)=0 wpp
#(b) s(1)=4 s(0)=s(2)=1 s(j)=0 wpp
#okresl nosniki z wykresu

x=-5:5;
y=zeros(size(x));
y(7)=1;
pp1=spline(x,y);
pp2=spline(x,[0,y,0]);
s=linspace(-5,5);
plot(s,ppval(pp1,s),";not-a-k;",s,ppval(pp2,s),";herm;",x,y,"+")
pause(1);

y(7)=4;y(6)=y(8)=1;
pp1=spline(x,y);
pp2=spline(x,[0,y,0]);
s=linspace(-5,5);
plot(s,ppval(pp1,s),";not-a-k;",s,ppval(pp2,s),";herm;",x,y,"+")
pause(1);


#zd 3 testy rzedu zbieznosci interp splajnami kubicznymi

function [ne,nel,ee,pp,x,s]=testerspc(f,a,b,N=20,df=[],M=1000)
#f - wskaznik do funkcji liczacej wartosci na wektorach 
#a,b- konce odcinkow 
#typ 0 - spljan not-a knot 
#N - ilosc wezlow siatki -1 
#M - rozmiar siatki do liczenia bledu dyskretnego 
#df - wskznik funkcji do pochodnej f (jak [] t osplajn not-a-knot
#Output
#ee blad na siatce M punktow rownoodlegluych na [a,b]
#ne,nel - norma sup ee i  norma L^2 ee (diskretna)
#pp -wspolczynniki splajnu
#x wezly interpolacji
#s -siatka M+1 punktowa rowwoodlegla na [a,b]
#przyklad
#ne=testerspc(@sin,0,pi)
#[ne,nel,ee,pp,x,s]=testerspc(@sin,0,pi,20,df=@cos,500);

x=linspace(a,b,N+1);
s=linspace(a,b,M+1);
y=f(x);
if(!isnull(df))
 pp=spline(x,y);
else
 #if(isvector(df))
 pp=spline(x,[df(a),y,df(b)]);
endif
%plot(s,ppval(pp1,s),s,ppval(pp2,s),x,y,"+")
ee=f(s)-ppval(pp,s); 
ne=norm(ee,'inf');
nel=norm(ee,2)*sqrt((b-a)/N); #sqrt( \sum_k ee(k)^2*h ) h=(b-a)/N

endfunction

S=10;
f=@(x) sin(x)
df=@(x) cos(x)
a=0;b=2*pi;
N=5;
ne=1;
for k=1:S,
 nep=ne;
 ne=testerspc(f,0,pi,N);
 printf("[%d] er=%g erp/er=",N,ne);
 if(k>1) printf("%g",nep/ne);endif
 printf("\n");
 N*=2;
endfor

f=@(x) 10./(1+x.*x);
a=0;b=2*pi;
N=5;
ne=1;
for k=1:S,
 nep=ne;
 ne=testerspc(f,0,pi,N);
 printf("[%d] er=%g erp/er=",N,ne);
 if(k>1) printf("%g",nep/ne);endif
 printf("\n");
 N*=2;
endfor



