echo on
########################################################################
#
#    Przykladowa sesja - przyklady zastosowan prostych narzedzi z algebry
#    liniowej oraz zastosowanie do znajdowania wielomianu najlepiej 
#    pasujacego do zadanych pktow - 'wprost' (rozwiazujemy
#    odpowiedni uklad rownan w sesie liniowego zad najmn. kwad.
#    i przy uzyciu funkcji polyfit - 
#    mozna tez jej uzyc do znalezienia wielomianu interp
#    Lagrange'a -  szczegolny przypadek - ilosc wezlow = st(wiel interp L)+1
#
#
########################################################################

#Tworzymy macierz
A=[ 1 2 ; 1 1 ]
f=[ 1 1]'
#Rozwiazujemy uklad rownan
x=A\f
#Sprawdzamy
norm(A*x-f,1)
#Uwarunkowanie macierzy?
cond(A)
#macierz odwrotna
Ai=inv(A)
#Sprawdzamy
norm(A*Ai-eye(size(A)),1)
norm(Ai*A-eye(size(A)),1)
#wartosci wlasne A
eig(A)
#i wektory wlasne
[V,D]=eig(A)
#Co gdy macierz prostokatna?
A=[ 1 2 3 ; 1 1 1]'
f= [ 1 1 1]'
#wtedy \ rozwiazuje uklad w sensie liniowych kwadratow
A \ f
#A gdy A nadokreslona?
A=A'
f=[ 3 2]'
#Oczywiste rozwiazanie to x=[ 1 1 0]'
x=A\f
#Sprawdzamy
norm(A*x-f,1)
#A roznica rozwiazan?
norm(A*(x- [1 1 0]'),1)
#czyli roznica w Ker(A)
#Czy da sie policzyc Ker(A)?
ker=null(A)
#sprawdamy
A*ker

#Co ze zle uwarunkowanymi macierzami?
#Macierz Hilberta
cond(hilb(10))
#Przetestujmy czy obliczona macierz odwrotna jest lepsza/gorsza od macierzy
#odwrotnej  obliczonej funkcja invhilb(N)
N=3;norm(invhilb(N)- inv(hilb(N)),1)
N=10;norm(invhilb(N)- inv(hilb(N)),1)
#Kiepsko!
#Sprawdzmy czy to sa w ogole odwrotnosci
#prawostronna
N=10;norm(hilb(N)*invhilb(N)- eye(N),1)
N=10;norm(hilb(N)*inv(hilb(N))- eye(N),1)
#Gorzej!
#Lewostronna
N=10;norm(invhilb(N)*hilb(N)- eye(N),1)
N=10;norm(inv(hilb(N))*hilb(N)- eye(N),1)
#Podobnie
#Czy lepiej rozwiazac uklad mnozac przez macierz odwrotna?
N=10;sol= 1:N;sol=sol';A=hilb(N);f=A*sol;
x=A\f;y=invhilb(N)*f;z=inv(hilb(N))*f;
norm(sol-x,1)
norm(sol-y,1)
norm(sol-z,1)
#Wielkiej roznicy nie ma

#Warto jeszcze znac svd, qr lu (tylko macierze kwadratowe)
#R, U dolnotrojkatne, U V Q - ortogonalne tzn Q'*Q=Id
# L-  dolnotrojkatna z jedynkami na diagonali
# S - diagonalna (dla prostokatnych tez)
A=[ 1 2 ; 3 4]
[U, S, V] = svd (A)
norm(A- U*S*V',1)
norm(U'*U-eye(size(U)),1)
norm(V'*V-eye(size(V)),1)
[Q,R]=qr(A)
norm(A-Q*R,1)
norm(Q'*Q-eye(size(Q)),1)
[L,U,P]=lu(A)
norm(P*A-L*U,1)

#Zastosujmy to co wiemy do zag znalezienia wielomianu najlepiej pasujacego
#do danej tabelki pktow - przypadek liniowy
#Sami: szukamy \min_{a,b} \sum|a*x(i) +b - y(i)|^2 = \sum|d*x(i) +c - y(i)|^2
x=1:10;y=1+x;
#Rozwiazanie a=[1 1]'
vander(x);
A=ans(:,9:10);
a=A\y'
#OK!
#A teraz przy pomocy funkcji polyfit
polyfit(x,y,1)
#to samo
#Funkcja tez moze sluzyc do interpolacji liniowej
polyfit(x(1:2),y(1:2),1)
#Czy wiecej wymiarowe (dlug danych = st wielomianu +1 np
polyfit(x,y,9)
#Tez OK, zauwazmy ze to wspolczynniki wielomianu st 9 - 
#dwa ostatnie sa rowne 1 a pozostale zera (no numeryczne zera)
