% Numerik D-MAVT, RW (Nipp), SS 06
% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HERE WE USE THE FUNCTIN jac DEFINED IN ANOTHER FILE %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Serie 3, A4
% b) Laplace equation on a square: check with exact solution
disp(['Aufgabe 4b']);
tol = 1e-5;
for k=[10 20 40]
    disp(['k=',num2str(k)]);
    A=model(k); b=rhs(k); exact = A\b;          % Needs funciton mode, rhs
    x0=ones(size(b)); % initial guess
    [x,its,cpu]=jac(A,b,x0,tol);                % Needs function jac
    disp(['Jacobi:         relative error = ',num2str(norm(exact-x)/norm(exact))]);
end
% MATLAB OUTPUT
%Aufgabe 4b
%k=10
%Jacobi:         relative error = 0.00015037
%k=20
%Jacobi:         relative error = 0.00071516
%k=40
%Jacobi:         relative error = 0.0030596



%c) Test solution of linear system using matrices A,B,C
disp(['Aufgabe 4c']);
tol = 1e-5;
A=[2 1 0;      % TEST WITH MATRIX A
    1 2 1;
    0 1 2];
r=ones(3,1);
exact=A\r;
x0 =ones(size(r));
[x,its,cpu]=jac(A,r,x0,tol);
disp(['Jacobi:            relative error = ',num2str(norm(exact-x)/norm(exact))]);
% MATLAB OUTPUT
%Aufgabe 4c
%Jacobi:            relative error = 3.3036e-06

tol =1e-5;
B=[1 -1 0;      % TEST WITH MATRIX B
    5 2 0;
    0 0 1];
r=ones(3,1);
exact=B\r;
x0=ones(size(B));
[x,its,cpu]=jac(B,r,x0,tol);
%MATLAB OUTPUT
% NO CONVERGENCE !!!!

tol =1e-5;
C=[3 2 1;      % TEST WITH MATRIX B
    2 3 2;
    1 2 3];
r=ones(3,1);
exact=C\r;
x0=ones(size(C));
[x,its,cpu]=jac(C,r,x0,tol);
%MATLAB OUTPUT
% NO CONVERGENCE !!!!




