
#zad 1-4 trywialne

#zad 5
#bledy dla 2 m grnotr U1=same 1 na diag i ponadiag
# U2=1 na diagonali ; -1 nad diagonala
format long e
nn=[4,16,23,52,60,100];

for k=1:length(nn),
 n=nn(k)
 U1=triu(ones(n,n));
 U2=-U1+2*eye(n,n);
 s=rand(n,1);ns=norm(s)
 f1=U1*s;nf1=norm(f1);
 f2=U2*s;nf2=norm(f2);
 x1=U1\f1;
 x2=U2\f2;
 e1=norm(x1-s)/ns
 r1=norm(U1*x1-f1)/nf1
 e2=norm(x2-s)/ns
 r2=norm(U2*x2-f2)/nf2
endfor


#zad 6 jak 5 ale dla m Hilberta Hn=(1/(i+j-1))_{i,j=1}^n
printf("testy dla m Hilberta\n");
format long e
nn=[4,16,23,30,52,60];

for k=1:length(nn),
 n=nn(k)
 H=hilb(n);
 s=rand(n,1);ns=norm(s);
 f1=H*s;nf1=norm(f1);
 x1=H\f1;
 e1=norm(x1-s)/ns
 r1=norm(H*x1-f1)/nf1
endfor

#zad 13
#Lu dla m 3diagonalnej bez permutacji wierszy

function x=solve3diag(a,b,c,f)
#Funkcja rozwiazuje ukald AX=f (o ile rozklad LU A istnieje)
#Wejscie:
#macierz A 3-diagonalna z wektorem a jako gl diagonala
#    b- poddiagonala a c naddiagonala
#f- wektor prawej strony
#Wyjscie:
#x- rozwiazanie Ax=f

n=length(a);
x=zeros(n,1);
#rozklad elementy l(k) w b a U w a (c sie nie zmienia)
for k=1:n-1,
 b(k)=b(k)/a(k);
 a(k+1)=a(k+1)-b(k)*c(k);
 f(k+1)=f(k+1)-b(k)*f(k);
endfor

x(n)=f(n)/a(n);
for k=(n-1):-1:1,
 x(k)=(f(k)-c(k)*x(k+1))/a(k);
endfor

endfunction

printf("\n testy dla A 3-diagonalnej \n");

nn=[10,100,1000,1000];

for k=1:length(nn),
 n=nn(k)
#a=2*ones(n,1);b=c=-ones(n-1,1);
a=rand(n,1);b=c=rand(n-1,1);
A=diag(a)+diag(b,-1)+diag(c,1);
s=rand(n,1);f=A*s;
tic;y=A\f;t1=toc
e1=norm(y-s)
r1=norm(A*y-f)
tic;x=solve3diag(a,b,c,f);t2=toc
e2=norm(x-s)
r2=norm(A*x-f)
endfor


