% ==========
% Apartado I
% ==========
%
%
% * * * EJERCICIO 1 * * *
A=[ 112,84,48,28; 84,81,77,-7; 48,77,146,-83; 28,-7,-83,85];
b=[-204;-186;-193;38];
x=A\b
eig(A)
disp('Factorizacion LU'),pause
[L,U1,P]=lu(A)
y=L\(P*b)
U1\y
disp('Factorizacion Choleski'),pause
U2=chol(A)
y=(U2')\b
U2\y
% ===========
% Apartado II
% ===========
%
%
% * * * EJEMPLO 1 * * *
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
KJ=D\(L+U) ,cJ=D\b 		% Método de Jacobi
KG=(D-L)\U ,cG=(D-L)\b	% Método de Gauss-Seidel
w=0.5 ;% Método de relajación con w=0.5
KS=(D-w*L)\((1-w)*D+w*U), cS=(D-w*L)\(w*b)
%
% Datos del ejercicio 2
%
%AA=3*eye(4);AA(1:3,2:4)=AA(1:3,2:4)+eye(3);AA(2:4,1:3)=AA(2:4,1:3)+eye(3)
AA=2*eye(4)+ones(4);AA=triu(AA,-1);AA=tril(AA,1)
bb=[-4;3;0;2];
D=diag(diag(AA));
L=-tril(AA,-1);
U=-triu(AA,1);
KJ=D\(L+U) ,cJ=D\bb 		% Método de Jacobi
KG=(D-L)\U ,cG=(D-L)\bb	% Método de Gauss-Seidel
w=0.5 ;% Método de relajación con w=0.5
KS=(D-w*L)\((1-w)*D+w*U), cS=(D-w*L)\(w*bb)

%
% * * * EJERCICIO 2 * * *
% ##################### Inicio del fichero ####################
function [K,c]=mitera(A,b,metodo,w)
%
%function [K,c]=mitera(A,b,metodo,w)
%Obtiene las matrices de iteración de los diferentes métodos
%
%Entrada
%	A,b         matrices
%	metodo      'jacobi','gauss','seidel','gauss-seidel','relajacion'
%   w           parametro de relajacion (si es necesario)
%Salida
%	K,c         Matrices del método iterativo x=Kx+c
%

% Verificacion de los argumentos de entrada
if nargin<3;error('Numero insuficiente de argumentos');end
[fa,ca]=size(A);if fa~=ca;error('La matriz no es cuadrada');end
[fb,cb]=size(A);if fa~=fb || cb~=1;error('El vector no es congruente');end
D=diag(diag(A));if any(diag(D)==0);error('matriz D no invertible');end
L=-tril(A,-1);
U=-triu(A,1);
switch lower(metodo)
    case {'jacobi'} 
        K=D\(L+U);c=D\b; 		% Método de Jacobi
    case {'gauss','seidel','gauss-seidel'}
        K=(D-L)\U;c=(D-L)\b;	% Método de Gauss-Seidel
    case {'relajacion'}
        if nargin<4;error('Numero insuficiente de argumentos');end
        K=(D-w*L)\((1-w)*D+w*U);c=(D-w*L)\(w*b);% Método de relajación
    otherwise
        error('método desconocido');
end
% ##################### fin del fichero ####################
%
% Ejecución de la función
%
[Kj,cj]=mitera(A,b,'Jacobi')
[Kg,cg]=mitera(A,b,'gauss')
[Ks,cs]=mitera(A,b,'relajacion',0.5)
%
% * * * EJERCICIO 3 * * *
w=0:0.1:2;rho=zeros(size(w));
N=length(w);
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
for i=1:N;
    K=(D-w(i)*L)\((1-w(i))*D+w(i)*U);rho(i)=max(abs(eig(K)));
end
plot(w,rho);title('Metodo de relajacion');grid on
xlabel('Parametro de relajacion');ylabel('Radio espectral')