% Numerik D-MAVT, RW (Nipp), SS06
% Serie 3, A1
format compact; 
% to check convergence look at spectral radius of iteration matrix T
A=[2 1 0; 1 2 1; 0 1 2]; 			% the matrices to test
B=[1 -1 0; 5 2 0; 0 0 1]
C=[3 2 1; 2 3 2; 1 2 3];
DA=diag(diag(A)); LA=tril(A)-DA; RA=triu(A)-DA;      % A=DA+LA+RA
DB=diag(diag(B)); LB=tril(B)-DB; RB=triu(B)-DB;      % B=DB+LB+RB
DC=diag(diag(C)); LC=tril(C)-DC; RC=triu(C)-DC;      % C=DC+LC+RC
TAj=-inv(DA)*(LA+RA); TAgs=-inv(DA+LA)*RA;	% x^{k+1}=T*x^k+c jacobi & gauss-seidel
TBj=-inv(DB)*(LB+RB); TBgs=-inv(DB+LB)*RB;	
TCj=-inv(DC)*(LC+RC); TCgs=-inv(DC+LC)*RC;	
% if rho < 1, we have convergence (assume (I-T) is regular)
rhoAj=max(abs(eig(TAj))), rhoAgs=max(abs(eig(TAgs))),
rhoBj=max(abs(eig(TBj))), rhoBgs=max(abs(eig(TBgs))),
rhoCj=max(abs(eig(TCj))), rhoCgs=max(abs(eig(TCgs))),
% output:
% rhoAj =
%    0.70710678118655
% rhoAgs =
%    0.50000000000000
% rhoBj =
%    1.58113883008419
% rhoBgs =
%    2.50000000000000
% rhoCj =
%    1.12409377442300
% rhoCgs =
%    0.60831218151474


% Numerik D-MAVT, RW (Nipp), SS06
% Serie 3, A2
A=[3 2 1; 2 3 2; 1 2 3];
D=diag(diag(A)); L=tril(A)-D; R=triu(A)-D;
% a) plot the spectral radius of T
ws=.01:.01:2-.01;% omegas
rho=[]; % vector to store spectral radii 
for i=1:length(ws)
  w=ws(i);
  T=inv(D+w*L)*(-w*R+(1-w)*D);	% iteration matrix
  rho=[rho max(abs(eig(T)))];
end
figure; 	% plot
plot(ws,rho); title('Spectral Radius for SOR'); xlabel('\omega'); 
ylabel('\rho(T(\omega))');
% b) 
w=0.4; T=inv(D+w*L)*(-w*R+(1-w)*D); nT=norm(T);
disp(['w=',num2str(w),': ||T||_2=',num2str(nT),', ',num2str(ceil(-10*log(10)/log(nT))),' iterations']);
w=1.4; T=inv(D+w*L)*(-w*R+(1-w)*D); nT=norm(T);
disp(['w=',num2str(w),': ||T||_2=',num2str(nT),', ',num2str(ceil(-10*log(10)/log(nT))),' iterations']);

% output:
% w=0.4: ||T||_2=0.90294, 226 iterations
% w=1.4: ||T||_2=1.1869, -134 iterations



function [x, iter, time] = jac(A,b, x0, tol)

% jacobi to solve linear system A*x=b
% in:	A	matrix of linear system
%	b	right hand side of system
%	x0	initial guess
%	tol	tolerance for convergence
% out:	x	approximative solution 
%	iter	number of iterations needed
% 	time	cpu time needed
% auth: fab@math.ethz.ch

tic;		% start measuring time
D=diag(diag(A)); L=tril(A)-D; R=triu(A)-D;	% decompose matrix
S=-(L+R);
x=x0;	% start with initial guess
maxIter=2e3; iter=0;
while (iter<maxIter)
  iter=iter+1;
  xold=x;
  x=D\(S*x+b);	% Jacobi iteration without inverse
  if (norm(x-xold)/norm(x)<tol) break; end	% converged!
end
time=toc;	% stop measuring time
if (iter==maxIter) error(['no convergence after ',int2str(iter),' iterations']); end;


function [x, iter, time] = sor(A,b, x0, tol,w)

% successive over relaxation to solve linear system A*x=b
% in:	A	matrix of linear system
%	b	right hand side of system
%	x0	initial guess
%	tol	tolerance for convergence
% 	w=1	relaxation parameter omega 
% out:	x	approximative solution 
%	iter	number of iterations needed
% 	time	cpu time needed
% auth: fab@math.ethz.ch

if (nargin==4) w=1; disp('Gauss-Seidel (omega=1)'); end

tic;		% start measuring time
D=diag(diag(A)); L=tril(A)-D; R=triu(A)-D;	% decompose matrix
Tl=D+w*L;
Tr=-w*R+(1-w)*D;
x=x0;	% start with initial guess
maxIter=2e3; iter=0;
while (iter<maxIter)
  iter=iter+1;
  xold=x;
  x = Tl\(Tr*xold+w*b);
  if (norm(x-xold)/norm(x)<tol) break; end	% converged!
end
time=toc;	% stop measuring time
if (iter==maxIter) error(['no convergence after ',int2str(iter),' iterations']); end;




% Numerik D-MAVT, RW (Nipp), SS06
% Serie 3, A3
% a) Laplace equation on a square: check with exact solution
n=10;
for tol=[1e-4 1e-5 1e-6 1e-7 1e-8]
  disp(['tol=',num2str(tol)]); 
  A=model(n); b=rhs(n); exact=A\b;
  x0=ones(size(b));	% initial guess
  [x,its,cpu]=jac(A,b,x0,tol);  
  disp(['Jacobi:       relative error = ',num2str(norm(exact-x)/norm(exact))]);
  [x,its,cpu]=sor(A,b,x0,tol);
    disp(['Gauss-Seidel: relative error = ',num2str(norm(exact-x)/norm(exact))]);
  for omega=[.4 1.4]
    [x,its,cpu]=sor(A,b,x0,tol,omega);
    disp(['SOR (w=',num2str(omega),'):  relative error = ',num2str(norm(exact-x)/norm(exact))]);
  end
end
% output:
% tol=0.0001
% Jacobi:       relative error = 0.0015021
% Gauss-Seidel (omega=1)
% Gauss-Seidel: relative error = 0.0006676
% SOR (w=0.4):  relative error = 0.0031985
% SOR (w=1.4):  relative error = 0.00018876
% tol=1e-05
% Jacobi:       relative error = 0.00015037
% Gauss-Seidel (omega=1)
% Gauss-Seidel: relative error = 7.1123e-05
% SOR (w=0.4):  relative error = 0.00032198
% SOR (w=1.4):  relative error = 2.2821e-05
% tol=1e-06
% Jacobi:       relative error = 1.5054e-05
% Gauss-Seidel (omega=1)
% Gauss-Seidel: relative error = 6.6907e-06
% SOR (w=0.4):  relative error = 3.2413e-05
% SOR (w=1.4):  relative error = 1.9402e-06
% tol=1e-07
% Jacobi:       relative error = 1.507e-06
% Gauss-Seidel (omega=1)
% Gauss-Seidel: relative error = 7.1278e-07
% SOR (w=0.4):  relative error = 3.1659e-06
% SOR (w=1.4):  relative error = 2.3457e-07
% tol=1e-08
% Jacobi:       relative error = 1.5087e-07
% Gauss-Seidel (omega=1)
% Gauss-Seidel: relative error = 6.7053e-08
% SOR (w=0.4):  relative error = 3.187e-07
% SOR (w=1.4):  relative error = 1.9942e-08

% b) Laplace equation on a square: compare different sizes
tol=1e-5;
for n=[10 20 40]
  disp(['n=',int2str(n)]); 
  A=model(n); b=rhs(n); 
  x0=ones(size(b));	% initial guess
  [x,its,cpu]=jac(A,b,x0,tol);  
  disp(['Jacobi:       ',int2str(its),' iterations, ',num2str(cpu),' seconds']);
  [x,its,cpu]=sor(A,b,x0,tol);
  disp(['Gauss-Seidel: ',int2str(its),' iterations, ',num2str(cpu),' seconds']);
  for omega=[.4 1.4]
    [x,its,cpu]=sor(A,b,x0,tol,omega);
    disp(['SOR (w=',num2str(omega),'):  ',int2str(its),...
          ' iterations, ',num2str(cpu),' seconds']);
  end
end

% output:
% n=10
% Jacobi:       137 iterations, 0.040064 seconds
% Gauss-Seidel (omega=1)
% Gauss-Seidel: 75 iterations, 0.021961 seconds 
% SOR (w=0.4):  257 iterations, 0.071199 seconds 
% SOR (w=1.4):  31 iterations, 0.010036 seconds 
% n=20
% Jacobi:       523 iterations, 6.5363 seconds 
% Gauss-Seidel (omega=1)
% Gauss-Seidel: 288 iterations, 3.9066 seconds 
% SOR (w=0.4):  949 iterations, 12.6415 seconds 
% SOR (w=1.4):  134 iterations, 1.8531 seconds 
% n=40
% Jacobi:       1777 iterations, 583.2934 seconds 
% Gauss-Seidel (omega=1)
% Gauss-Seidel: 996 iterations, 333.846 seconds 
% SOR (w=0.4):  3132 iterations, 1040.3408 seconds 
% SOR (w=1.4):  481 iterations, 160.9133 seconds 

% d) check if Jacobi and Gauss-Seidel converge for matrices of A1
b=ones(3,1); x0=b; tol=1e-10;
A=[2 1 0; 1 2 1; 0 1 2]; 
[x,its,cpu]=jac(A,b,x0,tol); 
disp(['A: Jacobi:      ',int2str(its),' iterations']);
[x,its,cpu]=sor(A,b,x0,tol); 
disp(['A: Gauss-Seidel: ',int2str(its),' iterations']);
B=[1 -1 0; 5 2 0; 0 0 1];
[x,its,cpu]=jac(B,b,x0,tol); 
disp(['B: Jacobi:      ',int2str(its),' iterations']);
[x,its,cpu]=sor(B,b,x0,tol); 
disp(['B: Gauss-Seidel: ',int2str(its),' iterations']);
C=[3 2 1; 2 3 2; 1 2 3];
[x,its,cpu]=jac(C,b,x0,tol); 
disp(['C: Jacobi:      ',int2str(its),' iterations']);
[x,its,cpu]=sor(C,b,x0,tol); 
disp(['C: Gauss-Seidel: ',int2str(its),' iterations']);
% output:
% A: Jacobi:      71 iterations
% A: Gauss-Seidel: 3 iterations
% ??? Error using ==> jac
% no convergence after 2000 iterations
% ??? Error using ==> sor
% no convergence after 2000 iterations
% ??? Error using ==> jac
% no convergence after 2000 iterations
% C: Gauss-Seidel: 47 iterations
% check with A2b
A=[3 2 1;2 3 2;1 2 3];
w=0.4; [x,its,cpu]=sor(A,b,x0,tol,w); 
disp(['SOR (',num2str(w),'): ',int2str(its),' iterations']);
w=1.4; [x,its,cpu]=sor(A,b,x0,tol,w); 
disp(['SOR (',num2str(w),'): ',int2str(its),' iterations']);
% output:
% SOR (0.4): 182 iterations
% SOR (1.4): 37 iterations
