
#zad 1 test \ jako solvera dla LZNK

A=[1 1; 1 3; 2 -1];
x=[1;2];
f=A*x;
y=A\f
norm(A*y-f)
norm(x-y)
#teraz dualne zadanie tzn LZNK z A^T
x1=[1;-1;2]
g=A'*x1;
y1=A'\g
norm(A'*y1-g)
norm(x1-y1)
#dlaczego? y ma min norme 2ga wsrod rozwiazan!
norm(y1)
norm(x1)

#zad 2 zad 1 ale za pomoca rozkladu qr (tylko regularne LZNK)
#A f x z zad 1
[Q,R]=qr(A)
norm(Q'*Q-eye(size(Q)),'fro')
norm(Q*Q'-eye(size(Q)),'fro')
norm(A-Q*R,'fro')
g=Q'*f;
R1=R(1:2,1:2);
w=R1\(g(1:2))
norm(x-w)
norm(A*w-f)
norm(w-y) #czy rozwiazania identyczne?

#f poza R(A)
f=[1;7;-10]
x=A\f; #przyjmijmy ze t odokladne rozw LZNK..
g=Q'*f;
y=R1\(g(1:2));
norm(A*x-f,2)
norm(A*y-f,2)
norm(R*y-g,2)
norm(x-y,2) 


#zad 3 uklad r normalnych f=[1;7;-10]
B=A'*A;
z=B\(A'*f)
norm(A*z-f)
norm(x-z)


#zad 4 \ a QR dla AX=f z A  m vanderemonde'a (kwadratowej nxn)

for n=10:5:30,
 x=rand(n,1)-0.5; #n wezlow rownoodleglych z (0,1)
 A=vander(x);#A m vandermonde'a dla wezlow x
 f=A*x; 
 tic;y=A\f;tb=toc;
 tic;[Q,R]=qr(A);w=R\(Q'*f);tlu=toc;
 tic;[L,U,P]=lu(A); z=U\(L\(P*f));tqr=toc;
 printf("\nOperator backslash[%d], t=%d er=%7.6e reser=%g \n",n,tb,norm(x-y),norm(A*y-f));
 printf("Rozklad QR        [%d],t=%d er=%10.6e reser=%g\n",n,tqr,norm(x-w),norm(A*w-f));
 printf("Rozklad LU        [%d],t=%d  er=%10.6e reser=%g\n",n,tlu,norm(x-z),norm(A*z-f));
endfor


#zad 5 czy H^2=I i czy H==H'
#i 6 wyzeruj skladowe poza pierwsza  losowy wektora  za pomocac m Householdera
#dodatkowo test tw z wykladu u,v losowe o tej samej normie 2-giej czy H u=v dla wektora Hous. v-u

function y=househ(x,w,knw=w'*w)
#funkcja mnozaca x przez mac householdera z wektorem H. w
#Input x wektor w -wektor housh. kknw(opcjonalnie) kwadra normy 2-giej w
#output 
#y=H_w*x 
 y=x-(2/knw)*w*(w'*x);
endfunction

n=10;
u=2*rand(n,1)-1;
nu=norm(u,2);
w=u; w(1)=w(1)-nu;
nw=w'*w;

y=househ(u,w,nw)
norm(y(2:n),2)

u=rand(n,1)-0.5;u=u/norm(u);
v=rand(n,1)-0.5;v=v/norm(v);
w=u-v;kwn=w'*w;
y=househ(u,w);
norm(y-v)
z=househ(v,w);
norm(z-u)

#wlasnosci H
u=rand(n,1)-0.5;w=rand(n,1)-0.5;
#czy ||Hu||_2=||u||_2?
norm(househ(u,w))-norm(u)

#czy H=H^{-1}

norm(u-househ(househ(u,w),w))

#czy H'=H
v=rand(n,1)-0.5;
v'*(househ(u,w))-(househ(v,w))'*u


#zad najlepiej pasujaca krzywa  ax^2+by^2-1=0
#do danych punktow

function a=krzywaps(x,y)
 n=length(x);
 f=ones(n,1);
 A=[x.*x,y.*y];
 a=A\f;
endfunction 

x=linspace(-.45,.45,20)';
y=sqrt(1-4*x.*x)+ 0.2*(rand(size(x))-0.5);#punkty blisko elipsy 4x^2+y^2=1
a=krzywaps(x,y)
printf("Elipsa %f*x^2+%f*y^*=1\n",a(1),a(2));
bl= a(1)*x.*x+a(2)*y.*y-1;
blad=sum(bl.*bl);
printf("Blad: sum_k |%3.1f*x^2+%3.1f*y^*-1|^2=%3.3f\n",a(1),a(2),blad);
s=linspace(-.5,.5);
plot(x,y,"+r",s,sqrt(1-4*s.*s),";4x^2+y^2=1;",s,sqrt(1-a(1)*s.*s)/a(2),";krzywa pasujaca;")


