


#zad 1 \ rozwiazuje LZNK
printf("Zad 1 backslash rozwiazuje LZNK\n");
A=[1 1;1 0;1,-1];
b=[1;1;1]
printf("czy wychodzi [1;0]?\n");
A\b #czy wersor?

#zad 2 LZNK na 3 sposoby A z z1
printf("Zad 2 - 3 metody  rozwiazywania LZNK z z1\n");
b=[1;2;2];
##\
x=A\b
#ukl r normalnych
B=A'*A;g=A'*b;
y=B\g
printf("czy x=yq?\n");
norm(x-y)
#OK

#poprzez qr na piechote
[m,n]=size(A);
[Q,R]=qr(A);
g1=Q'*b;
g=g1(1:n);
R=R(1:n,:);
z=R\g
printf("czy x=z?\n");
norm(x-z)
#blad \|Ax-b|_2
er=norm(A*x-b)
er1=norm(g1(n+1:m))
#czy rowne?
printf("czy bledy rowne sobie?\n");
abs(er-er1)

#poprzez gramma Schmidta 
#zamieszczam funkcje...


function [Q,R]=gramsch(A)
#ortogonalizacja Gramma Schmidta dla macierzy N kolumnowej
#Input 
#A - macier m x n (m>=n)
#Output 
#Q macierz o kolumnach ortonormalnych m x n
#R macierz nieosobliwa gornotrojkatna n x n

tolfl=1e2*eps; #tolerance for fl if R(k,k)<tolfl -program breaks
[m,n]=size(A);
if(m<n) 
 printf("\n No of columns > no of rows \n");
 return;
endif
 Q=zeros(m,n);
 R=zeros(n,n);
 for k=1:n,
  q=A(:,k);
  R(1:k-1,k)=Q(:,1:k-1)'*q;
  q=q - Q(:,1:k-1)*R(1:k-1,k);
  R(k,k)=sqrt(q'*q);
  if(R(k,k)>tolfl) 
    q=q/R(k,k);
  else 
    printf("\n Warning: Numerical singularity ||q(%d)||_2=%g<%g!\n",k,R(k,k),tolfl);
    return;
  endif
  Q(:,k)=q;
 end
endfunction

[Q1,R1]=gramsch(A);
#testujemy
norm(A-Q1*R1,1)
norm(Q1'*Q1-eye(size(R1)),1)
w=R1\(Q1'*b);
printf("czy x=w?\n");
norm(x-w)
#blad \|Ax-b|_2
er2=norm(A*w-b)
#czy rowne?
printf("czy bledy rowne sobie?\n");
abs(er-er2)
#w powinno byc rowne z?
norm(w-z)
#porownajmy Q i Q1


#zad 3 punkty najlepiej pasujace do krzywej ax^2+by^2=1
printf("Zad 3 -punkty najlepiej pasujace do krzywej ax^2+by^2=1\n");
n=20;

x=linspace(-1,1,n)';
y=sqrt(1-x.*x);
A=[x.*x,y.*y];f=ones(n,1);
printf("czy z=[1;1]?\n");
z=A\f

#teraz z zaburzeniem
y=sqrt(1-x.*x)+(rand(n,1)-0.5)*1e-3;
A(:,2)=y.*y;
printf("czy z bliski [1;1]?\n");
z=A\f

s=linspace(-1,1);
w=sqrt((1-z(1)*s.*s)/z(2));
plot(x,y,"+",s,w)
pause(1);
#zwiekszmy zaburzenie
y=sqrt(1-x.*x)+(rand(n,1)-0.5)*1e-1;
A(:,2)=y.*y;

z=A\f

s=linspace(-1,1);
w=sqrt((1-z(1)*s.*s)/z(2));
plot(x,y,"+",s,w)
pause(1);

#z4 p Householdera
printf("Zad 4 - macierz householdera wlasnosci\n");

function y=H(x,w,knw=w'*w)
#Funkcja mnozaca wektor/macierz przez H macierz Householdera 
#z wektorem Householdera w
#
#Input
#w wektor householdera
#x - wektor (macierz)
#knw - kwadrat normy w (opcja)
#Output
#y- Hx (x przemnozony przez macierz householdera z wektorem H. w
y=x-(2/knw)*w*(w'*x);
endfunction

#testujemy normy czy mnozenie wektora/macierzy przez mac. H Houdsholdera 
#zmienia norme 2 czy Frobeniusa
printf("czy mnozenie prez m Householdera zmienia normy 2 i Frobeniusa\n");
n=20;
printf("Powinny byc same zera\n");
for k=1:30,
 x=rand(n,1);
 w=rand(n,1);
 A=rand(n,n);
 norm(x,2)-norm(H(x,w),2) 
 norm(A,2)-norm(H(A,w),2) 
 norm(A,2)-norm(H(A',w)',2) 
 norm(A,'fro')-norm(H(A,w),'fro')  
endfor

#testujemy czy Hx=y dla x,y o tych samych normach 2 z wek. Househ. w=x-y
printf("czy dla x,y o tej samej normie 2 i dla  H z wektorem H. w=x-y spelnia y=Hx?\n");
n=10;
x=rand(n,1);
y=zeros(n,1);
y(1)=norm(x,2);
printf("Powinny wyjsc zera\n");
norm(y-H(x,x-y))
norm(-y-H(x,x+y))

#zerowanie 1-szej kolumny - 2 mozliwosci
printf("zerowanie pierwszej kolumny macierzy A (elementy pod diagonala)\n");
n=5;
A=rand(n,n);
y=zeros(n,1);
y(1)=norm(A(:,1),2);
printf("Czy w pierwszej kolumnie są zera pod diagonala?\n");
H(A,A(:,1)+y)
H(A,A(:,1)-y)
