#!/usr/bin/octave

################################################################################
#zad 4 i 5
# zbieznosc interpolacji splajnowej dla podwajanej ilosci wezlow
#dla funkcji sin(x) na [-20,2*pi]
################################################################################
K=10;
a=-20;b=2*pi;
N=5;h=(b-a)/N;
z=linspace(a,b,10000);
w=linspace(a,b,1000);
e=ep=eh=ehp=0;
f=@(x) sin(x);
df=@(x) cos(x);
ff=f(z);
for k=1:K,
x=linspace(a,b,N);
y=f(x);
#spline interpolacyjny not-a-knot
pp=spline(x,y);
#spline hermitowski
ph=spline(x,[df(a),y,df(b)]);
#ehp=eh;
ep=e;ehp=eh;
#bledy miedzy splajnami a funkcja w normie sup (dyskretnej)
e=max(abs(ff-ppval(pp,z)));
eh=max(abs(ff-ppval(ph,z)));
#drukujemy bledy w normie supremum
printf("N=%d h=%g  e=%g ep/e=%g e/(h^4)=%g ",N,h,e,ep/e,e/(h^4));
printf("eh=%g ehp/eh=%g eh/(h^4)=%g \n",eh,ehp/eh,eh/(h^4));
#dla malej ilosci bledow rysujemy wykresy splajnow i funkcji
if (k<4) 
  plot(w,ppval(pp,w),";not-a-knot;",w,ppval(ph,w),";herm;",w,f(w),";f;",x,y,"*;wezly;");
  pause(3);
end
#A tu rysujmey wykresy bledow
#plot(w,ppval(pp,w)- f(w),";blad;");pause(1);
N*=2;
h/=2;
endfor
pause(3);

