
disp("testy metody potegowej");



function [x,r,it]=power17(A,tol=1e-12,x0=randn(columns(A),1))
%metoda potegowa dla macierzy symetrycznej
%x(k+1)=Ax(k)/||A*x(k)||_2
%warunek stopu kolejne ilorazy Rayleigha roznia sie o mniej niz zadana tolerancja TOL
%INPUT 
%A diagonalizowalna (symetryczna)
%tol - tolerancja - domyslnie 10^(-12)
%x0 - domyslnie losowy
%OUTPUT
%x - wektor wlasny dla w.w. (przyblizenie) 
%r - wartosci wlasnej (przyblizenie)
%Przyklad: 
%l1=20;M=10;[Q,R]=qr(randn(M+1,M+1)); A=Q*diag([l1,1:M])*Q';[x,r,it]=power17(A)
%[r-l1,norm(x-Q(:,1)),norm(x+Q(:,1))] %blad przyblizenia wartoci i wektora wlasnego 


x=x0/norm(x0);
Ax=A*x;
it=0;
r=x'*Ax;
x=Ax/norm(Ax);;
it++;
Ax=A*x;
rp=r;
r=x'*Ax;
x=Ax/norm(Ax);;
Reer=abs(r-rp);
 while(Reer>tol)
 it++;
 Ax=A*x;
 rp=r;
 r=x'*Ax;
 x=Ax/norm(Ax);
 Reer=abs(r-rp);
 end

end

function [x,r,it]=ipower17(A,a=0,tol=1e-12,x0=randn(columns(A),1))
%metoda potegowa odwrotna dla macierzy  symetrycznej z parametrem 
%x(k+1)=Ax(k)/||A*x(k)||_2
%warunek stopu kolejne ilorazy Rayleigha roznia sie o mniej niz zadana tolerancja TOL
%INPUT 
%A diagonalizowalna (symetryczna)
%tol - tolerancja - domyslnie 10^(-12)
%x0 - domyslnie losowy
%OUTPUT
%x - wektor wlasny dla w.w. (przyblizenie) 
%r - wartosci wlasnej (przyblizenie)
%Przyklad: 
%M=11;[Q,R]=qr(randn(M+1,M+1)); A=Q*diag([l1,1:M])*Q';[x,r,it]=ipower17(A)
%[r-1,norm(x-Q(:,2)),norm(x+Q(:,2))] %blad przyblizenia wartoci i wektora wlasnego 
%M=11;[Q,R]=qr(randn(M+1,M+1)); A=Q*diag([l1,1:M])*Q';[x,r,it]=ipower17(A,7.01)
%[r-7,norm(x-Q(:,8)),norm(x+Q(:,8))] %blad przyblizenia wartoci i wektora wlasnego  




A=A-a*eye(size(A));
x=x0/norm(x0);
Ax=A\x;
it=0;
r=x'*Ax;r=1/r+a;
x=Ax/norm(Ax);;
it++;
Ax=A\x;
rp=r;
r=x'*Ax;r=1/r+a;
x=Ax/norm(Ax);;
Reer=abs(r-rp);
 while(Reer>tol)
 it++;
 Ax=A\x;
 rp=r;
 r=x'*Ax;r=1/r+a;%x'*Ax to iloraz Raylegha dla (A-aI)^{-1}
 x=Ax/norm(Ax);
 Reer=abs(r-rp);
 end

end
echo("on");
%tworzymy losowa macierz ortogonalna M+1 x M+1
l1=20;M=10;[Q,R]=qr(randn(M+1,M+1)); 
%teraz macierz o wartosciach wlasnych l1,1,2,3,..,M
A=Q*diag([l1,1:M])*Q';
%obliczamy przyblizenia pary wlasnej metoda potegowa
[x,r,it]=power17(A)
%i testujemy blad
[r-l1,norm(x-Q(:,1)),norm(x+Q(:,1))] %blad przyblizenia wartoci i wektora wlasnego 

%powtarzamy to samo dla m odwrotnej potegowej z parametrem a=0
[x,r,it]=ipower17(A)
%testujemy czy iloraz Rayleigha zbiegl do najmniejszej w walsnej tzn 1
[r-1,norm(x-Q(:,2)),norm(x+Q(:,2))] %blad przyblizenia wartoci i wektora wlasnego 
%a teraz a bliska wartosci wlasnej 7
[x,r,it]=ipower17(A,7.01)
%patrzymy na blad
[r-7,norm(x-Q(:,8)),norm(x+Q(:,8))] %blad przyblizenia wartoci i wektora wlasnego  

%zawsze blad przyblizenia ww kilka rzedow wielksoci mniejszy niz norma bledy przyblizenia wektora wlasnego 

%wezmy macierz diagonalizowalna ale niesymetryczna

%tworzymy losowa macierz  M+1 x M+1
l1=20;M=10;C=randn(M+1,M+1); 
%teraz macierz o wartosciach wlasnych l1,1,2,3,..,M
A=C*diag([l1,1:M])*inv(C);
%obliczamy przyblizenia pary wlasnej metoda potegowa
[x,r,it]=power17(A)
%i testujemy blad
v=C(:,1)/norm(C(:,1));
[r-l1,norm(x-v),norm(x+v)] %blad przyblizenia wartoci i wektora wlasnego 
%ten sam poziom bledy dla wartoci i wektora wlasnego
%powtorz dla odwrotnej potegowej biorac np a=7.01 - porownaj wynik  z para wlasna

%zadanie 3
%testujemy czy biorac x0 wektor startowy ortogonalny do wektora walsnego dla dominujacj ww
%dostaniemy zbieznosc do wektora wlasnego dla 2giej co do wielkosci ww
%tworzymy losowa macierz ortogonalna M+1 x M+1
l1=123;l2=15;M=10;[Q,R]=qr(randn(M+1,M+1)); 
%teraz macierz o wartosciach wlasnych l1,1,2,3,..,M
A=Q*diag([l1,l2,2:M])*Q';
%obliczamy przyblizenia pary wlasnej metoda potegowa
x0=randn(M+1,1); 
%ortogonalizujemy wzgledem Q(:,1) wektora walsnego dla l1 
x=x0-(x0'*Q(:,1))*Q(:,1);
x'*Q(:,1) %czy zero?
MM=25;
rr=zeros(MM,5);
echo("off");
for k=1:14,
   Ax=A*x;
   r=Ax'*x;
   x=Ax/norm(Ax);
   rr(k,:)=[k,r,r-l1,r-l2,x'*Q(:,1)];
end
echo("on");
printf("k,r-l1,r-l2,x'*q1,\n");
rr(1:14,:)
%wyglada ze Ok zbiega do l2 2giej co do wielkosci w module ww
%iterujmy dalej
echo("off");
for k=15:MM,
   Ax=A*x;
   r=Ax'*x;
   x=Ax/norm(Ax);
   rr(k,:)=[k,r,r-l1,r-l2,x'*Q(:,1)];
end
echo("on");
printf("k,r-l1,r-l2,x'*q1,\n");
rr(15:MM,:)
disp("Niestety z powodu bledu zaookraglen x stracil ortogonalnosc od q1 i zbiegl znow do l1");
r-l1
eig(A)


%czy cos mozna zrobic? reortogonalizujemy co np 5 iteracji
x=x0-(x0'*Q(:,1))*Q(:,1);
x'*Q(:,1) %czy zero?
MM=40;
rr=zeros(MM,5);
echo("off");
for k=1:MM,
   Ax=A*x;
   r=Ax'*x;
   x=Ax/norm(Ax);
   if(mod(k,5)==0)
     %printf("k=%d reortogonalizumy x\n",k);
     x=x-(x'*Q(:,1))*Q(:,1);x=x/norm(x);
   end  
   rr(k,:)=[k,r,r-l1,r-l2,x'*Q(:,1)];
end
echo("on");
printf("k,r-l1,r-l2,x'*q1,\n");
rr
r-l2
%zbieglo do  2giej ww co do wielkosci (w module)
%moze starczy rzadziej ortogonalizowac? 
%prosze przetestowac z reortogonalizacja co np 10 lub 20 iteracji
%to zalezy od stosunku 2giej do 1szej ww (uszeregowanych wg wielkosci w module)
echo("off");