echo on

##########################################################################
# skrypt rozwiazujacy problem znaleznienia wielomianu 
# najlepszej aproksymacji w L^2(0,1) w bazie potegowej 
# uzywajacy funkcji L2approx
# rozwiazujac uklad normalny (z macierza Grama w bazie potegowej)
# najpierw definiujemy funkcje ktora znajdzie nam ten wielomian
# czyli policzy calki f(k+!)\int_0^1f(X)x^kdx k=0,..n  oraz potem 
# rozwi±ze uklad z macierza Grama - szczesliwie ze to akurat macierz Hilberta
# (oczywiscie w praktyce  powinien to byc osobny m-plik z 2 funkcjami - 
# druga subfunction)
############################################################################3

function p=L2app(nazwa_f,n)
#Funkcja wyznacza wspolczynniki wiel najlepszej aproksymacji 
# dla funkcji nazwa_f na odcinku [0,1] w normie L^2(0,1)
# nazwa_f - nazwa fnkcji potaci 
#function y=nazwa_f(x)

global L2appnazwa_f;
global L2stjedn;
L2appnazwa_f = nazwa_f;
L2stjedn=n;


      
f=zeros(n+1,1);

#Liczymy calki \int_0^1f(x)x^kdx
for k=0:n,
  L2stjedn=k;
  f(k+1)=quad("fxk",0,1);
end  

#tworzymy macierz Gramma czyli tu Hilberta
H=hilb(n+1);

#Rozwiazujemy uklad rownan
p=H\f;

#Musimy zmienic kolejnosc aby uzywac funkcji obslugujacych wielomiany
#w Octavie np polyval
p=p(n+1:-1:1);
endfunction


#def pomocniczej funkcji subfunkcji
function y=fxk(argx)
      global L2stjedn;
      global L2appnazwa_f;
      k=L2stjedn;
      a=1;x=argx;
      for s=1:k, a=x*a; endfor;
      y=feval(L2appnazwa_f,x);
       y=a*y;
endfunction



disp("L2 approximation - test");


#definiujemy funkcje testowa
function y=L2apptest(x)
 y = 1+x.*x;;
endfunction

#Wyznaczamy wspolczynnki wielomianu
p=L2app("L2apptest",4)

#OK ten sam wielomian bo f wiel st <=4! 0x^4+0x^3+ x^2 + 1 

#teraz dla sin(2*pi*x)
function y=L2apptest(x)
 y = sin(2*pi*x);
endfunction

p=L2app("L2apptest",4)

#Policzmy norme

#norma L2 (dyskretna) przybliza ciagla

disp('Norma pomiedzy funkcja a funkcja przybl');
mesh=linspace(0,1);
norm(polyval(p,mesh) - L2apptest(mesh),2)/length(mesh)
#OK! (dyskretna norma 2 aproksymuje norme L^2)

#A teraz dla funkcji indykatorowej [0,0.5]
function y=L2apptest(x)
 y = (x<0.5);
endfunction

p=L2app("L2apptest",6);


mesh=linspace(0,1);
norm(polyval(p,mesh) - L2apptest(mesh),2)/length(mesh)
#Zupelnie Ok!

echo off
#Wykres?
plot(mesh,L2apptest(mesh),mesh,polyval(p,mesh));

sleep(3);

#A dla wiekszego stopnia?
p=L2app("L2apptest",11);
plot(mesh,L2apptest(mesh),mesh,polyval(p,mesh));
#Skad ten warning i co oznacza?
