

%lab2.m
graphics_toolkit("gnuplot")

%rozwiazywanie ukladow r liniowych Ax=b
A=[1 1;-1 1];b=[2;0];
x=A\b
norm(A*x-b)
norm(x-[1;1])

A=[1 1;1 0.99];b=[2;1];
x1=A\b
A(2,2)=0.98;%zaburzamu o 1e-2 element na pozycji (2,2)
x2=A\b
norm(x1-x2)

%sprawdzamy rzklady
A=randn(4);A=A'*A +20*eye(size(A));
[l,u,p]=lu(A);
norm(p*A-l*u,1)
norm(l-tril(l),1)
norm(u-triu(u),1)
%czy U=D*L' dla m D diagonalnej z diag U
D=diag(diag(u));
norm(u-D*l',1)

R=chol(A);
norm(R'*R-A,1)
%jak R ma sie do L i D z PA=LU
sR=sign(diag(R));
sR=diag(sR);
norm(R'*sR-l*sqrt(D),1)

[Q,R]=qr(A);
norm(Q'*Q-eye(size(Q)),1)
norm(Q*R-A,1)


%rozwiaz 4 metodami backslash lu qr chol 
%uklad Ax=b dla A=A'>0 losowej roznych wymiarow

n=10;
for k=1:3,
  A=randn(n);A=A'*A;
  sol=randn(n,1);b=A*sol;
  x1=A\b;
  [l,u,p]=lu(A);x2=u\(l\(p*b));
  R=chol(A);x3=R\(R'\b);
  [Q,R]=qr(A);x4=R\(Q'*b);
  [norm(x1-sol),norm(x2-sol),norm(x3-sol),norm(x4-sol)]
   n=n*10;
end  

%porownajmy czas tych 4 algorytmow

n=10;
for k=1:3,
  A=randn(n);A=A'*A;
  sol=randn(n,1);b=A*sol;
  tic;x1=A\b;t(1)=toc;
  tic;[l,u,p]=lu(A);x2=u\(l\(p*b));t(2)=toc;
  tic;R=chol(A);x3=R\(R'\b);t(3)=toc;
  tic;[Q,R]=qr(A);x4=R\(Q'*b);t(4)=toc;
  [n,t]
  n=n*10;
  disp("1:wmiar;2 backslash, 3 lu; 4 chol; 5 qr");
end  


n=20;
M=8;
er=zeros(M,5);
for k=1:M,
  U1=triu(ones(n));
  U2=U1-2*eye(n);
  x=randn(n,1);x=x/norm(x);
  b1=U1*x;b2=U2*x;
  x1=U1\b1;x2=U2\b2;
  er(k,:)=[n,norm(x1-x),norm(b1-U1*x1),norm(x2-x),norm(b2-U2*x2)];
  n+=5;
endfor
 er
disp("W 4 tej kolumnie blad ||x-x2||_2 niepokojaco rosnie. Dlaczego?");
disp("z kolei blad residualny ||U2x2-b2||_2 na poziomie 1e-16...");
[cond(U1),cond(U2)] %uwarunkowanie U2 wielkie
