 

#zad 1 \ czy dziala

A1=[1 1;1 -1];
f=A1(:,1);
x=A1\f
norm(A1*x-f)
#A osobliwa
A=[1 1; 2 2]
f=A(:,1);
x=A\f
norm(A*x-f)


inv(A)*A

inv(A1)*A1

#macierz odwrotna

#czas dla losowej A
N=10
for k=1:4,
 N*=3
 A=randn(N,N);
 x=randn(N,1);
 f=A*x;iA=inv(A);
 tic;y=A\f;toc
norm(x-y,'inf')
 norm(A*y-f,'inf')
 tic;y=iA*f;toc
norm(x-y,'inf')
 norm(A*y-f,'inf')
 tic;y=inv(A)*f;toc
end


#macierz 2diag i 3diagonalne
N=6;

#macierz z 1 na diagonali i -1 na naddiagonali
U=diag(ones(N,1)) -diag(ones(N-1,1),1)
#odwrotna czy tez pasmowa? np 2 diagonalna?
iU=inv(U)
#A 3diagonalna
A=U+U'
#odwrotna czy pasmowa?
iA=inv(A)

#porownujemy czas rozwiazanai ukl r liniowych
#z U operatorem \ i przez mnozenie przez macierz odwrotna
#utworzona wczesniej explicite
for k=1:6,
 N*=3
 U=diag(ones(N,1)) -diag(ones(N-1,1),1);
 x=randn(N,1);
 f=U*x;iU=inv(U);
 tic;y=U\f;toc
 norm(x-y,'inf')
 norm(U*y-f,'inf')
 tic;y=iU*f;toc
 norm(x-y,'inf')
 norm(U*y-f,'inf')
 tic;y=inv(U)*f;toc
end

#3 diagonalna A=U+U'

#to samo dla 3diagonalnej
N=6
for k=1:5,
 N*=3
 U=diag(ones(N,1)) -diag(ones(N-1,1),1);
 A=U+U';
 x=randn(N,1);
 f=A*x;iA=inv(A);
 tic;y=A\f;toc
norm(x-y,'inf')
 norm(A*y-f,'inf')
 tic;y=iA*f;toc
norm(x-y,'inf')
 norm(A*y-f,'inf')
 tic;y=inv(A)*f;toc
end


#macierz Hilberta
M=30; 
condA=zeros(M,1);
nerr=nres=zeros(M,2);
for N=0:M,
 N++
 A=hilb(N);
 condA(N)=cond(A);
 x=randn(N,1);
 f=A*x;iA=inv(A);
 y=A\f;
 disp("bledy przy rozw backslashem");
 ner=nerr(N,1)=norm(x-y,'inf')
 nr =nres(N,1)=norm(A*y-f,'inf')
 y=iA*f;
 disp("bledy przy uzyciu inv(A)*f");
 ner=nerr(N,2)=norm(x-y,'inf')
 nr =nres(N,2)=norm(A*y-f,'inf')
end
dnerr=nerr;dnres=nres;
#by porownac...

#powtarzamy w single

#macierz Hilberta
#rozwiazujemy ukladz macierza Hilbera przez \ i nozenie przez m odwrotna
#liczymy blad i blad residualny tzn normy x-y i  Ay-f dla y - obliczonego rozwiazania
#x znane losowe rozwiazanie 
M=30; 
condA=zeros(M,1);
nerr=nres=zeros(M,2);
for N=0:M,
 N++
 A=single(hilb(N));
 condA(N)=cond(A);
 x=single(randn(N,1));
 f=A*x;iA=inv(A);
 y=A\f;
 disp("bledy przy rozw backslashem");
 ner=nerr(N,1)=norm(x-y,'inf')
 nr =nres(N,1)=norm(A*y-f,'inf')
 y=iA*f;
 disp("bledy przy uzyciu inv(A)*f");
 ner=nerr(N,2)=norm(x-y,'inf')
 nr =nres(N,2)=norm(A*y-f,'inf')
end
snerr=nerr;snres=nres;

[dnerr(1:7,:), snerr(1:7,:)]


