



N=10;
function [x,nr,it]=richardson(A,b,tau,x0=zeros(size(b)),tol=1e-6)
%metoda Richardsona x(n+1)=x(n)+tau*r(n) r(n)=b-A*x(n)
x=x0;
r=b-A*x;
nres=norm(r);

tol=tol*nres; 
maxit=100;
nr=zeros(maxit,2);
nr(1,:)=[0,nres];
it=0;
while((it<maxit)&&(nres>tol))
 x=x+tau*r;
 r=b-A*x;
 nres=norm(r);
 it++;
 nr(it,:)=[it,nres];
end
 nr=nr(1:it,:);
endfunction




N=10;
Q=rand(N,N);
[Q,R]=qr(Q);
clear R;
D=diag(N:-1:1);
A=Q*D*Q';
eig(A)
s=rand(N,1);
b=A*s;


tau=1/6; %< 2/10
[x,nr,it]=richardson(A,b,tau)

tauo=2/(1+10);

[x,nres,it]=richardson(A,b,tauo)



function [x,nr,it]=jakobi(A,b,x0=zeros(size(b)),tol=1e-6)
%metoda Jakobi Dx(n+1)+(L+U)x(n)=b
%              Dx(n+1)+(A-D)x(n)=b czyli x(n+1)=x(n) +D^{-1}(b-Ax(n))  
%x(n+1)=x(n)+D^{-1}*tau*r(n) r(n)=b-A*x(n)
x=x0;
r=b-A*x;
nres=norm(r);
d=diag(A);di=1./d;
tol=tol*nres; 
maxit=100;
nr=zeros(maxit,2);
nr(1,:)=[0,nres];
it=0;
while((it<maxit)&&(nres>tol))
 x=x+di.*r;
 r=b-A*x;
 nres=norm(r);
 it++;
 nr(it,:)=[it,nres];
end
 nr=nr(1:it,:);
endfunction

function [x,nr,it]=gs(A,b,x0=zeros(size(b)),tol=1e-6)
%metoda Gauss-Seidel (D+L)x(n+1)+U x(n)=b
%              (D+L)x(n+1)+(A-(D+L))x(n)=b czyli x(n+1)=x(n) +(D+L)^{-1}(b-Ax(n))  
%x(n+1)=x(n)+(D+L)^{-1}*tau*r(n) r(n)=b-A*x(n)
x=x0;
r=b-A*x;
nres=norm(r);
DL=tril(A);
tol=tol*nres; 
maxit=100;
nr=zeros(maxit,2);
nr(1,:)=[0,nres];
it=0;
while((it<maxit)&&(nres>tol))
 x=x+DL\r;
 r=b-A*x;
 nres=norm(r);
 it++;
 nr(it,:)=[it,nres];
end

 nr=nr(1:it,:);
endfunction

N=10;
Q=rand(N,N);
[Q,R]=qr(Q);
clear R;
D=diag(N:-1:1);
A=Q*D*Q';
eig(A)
s=rand(N,1);
b=A*s;



tau=1/6; %< 2/10
[x,nr,it]=richardson(A,b,tau)

tauo=2/(1+10);

[x,nres,it]=richardson(A,b,tauo)
[x,nr,it]=gs(A,b)
[x,nr,it]=jakobi(A,b)




N=10;
Q=rand(N,N);
[Q,R]=qr(Q);
clear R;
D=diag(N:-1:1);
A=Q*D*Q';
eig(A)
s=rand(N,1);
b=A*s;



tau=1/6; %< 2/10
[x,nr,it]=richardson(A,b,tau)

tauo=2/(1+10);

[x,nres,it]=richardson(A,b,tauo)
[x,nr,it]=gs(A,b)
[x,nr,it]=jakobi(A,b)

%metoda potegowa

function [r,x,it]=potega(A,x0=rand(rows(A),1))
%metoda potegowa x=A*x/||Ax||_2
%INPUT
%A symetryczna
%x0 -wektor startowy  default losowy wektor
%OUTPUT
%iloraz Rayleigha r=x*A*x, 
%x przyblizenie wektora wlasnego dla najwiekszej w module ww
%it ilosc iteracji 

  maxit=100;
  tol=1e-14;
  x=x0/norm(x0);
  Ax=A*x; %A*x0
  rp=x'*Ax; %r0
  x=Ax/norm(Ax); %x1
  Ax=A*x; %A*x1
  r=x'*Ax; %r1
  it=1;
  while((it<maxit)&&(abs(r-rp)>tol))
   x=Ax/norm(Ax);
   Ax=A*x;
   rp=r;r=x'*Ax;
   it++;
  end
endfunction
B=[2,1;1,2]
[r,x,it]=potega(B)

[r,x,it]=potega(A)
D=diag((N:-1:1));
D(1,1)=100;
A=Q*D*Q';
[r,x,it]=potega(A)
[r-100,norm(x-Q(:,1)),norm(x+Q(:,1))]



function [r,x,it]=odwpot(A,a=0,x0=rand(rows(A),1))
%odwrotna metoda potegowa x=B*x/||Bx||_2
%dla B=(A-a*I)^(-1)
%INPUT
%A symetryczna
%a parametr
%x0 -wektor startowy  default losowy wektor
%OUTPUT
%iloraz Rayleigha r=x*A*x, 
%x przyblizenie wektora wlasnego dla najwiekszej w module ww
%it ilosc iteracji 
  A=A-a*eye(size(A));
  maxit=100;
  tol=1e-14;
  x=x0/norm(x0);
  iAx=A\x; %A\x0
  rp=x'*iAx; %r0
  x=iAx/norm(iAx); %x1
  iAx=A\x; %A*x1
  r=x'*iAx; %r1 
  it=1;
  while((it<maxit)&&(abs(r-rp)>tol))
   x=iAx/norm(iAx);
   iAx=A\x;
   rp=r;r=x'*iAx;
   it++;
  end
  r=1/r+a; %obliczamy przyblizenie ww A (bo r przyblizal ww (A-a)^(-1))
endfunction

[r,x,it]=odwpot(A)
[r,x,it]=odwpot(A,2.1)
[r,x,it]=odwpot(A,2.0001)

