% ============
% Apartado 1.1
% ============
%
%
% ##################### Inicio del fichero ####################
function R=opera(A,operacion,B)
%
%function function R=opera(A,operacion,B)
%Realiza operaciones con una y otra matriz dependiendo de la operacion
%Ejemplo de uso de la instrucción case
%
%Entrada
%	A,B         matrices
%	operacion   suma('+'),resta('-'),multiplicación('*') o inversión(#')
%Salida
%	R           resultado
%
% Verificacion de los argumentos de entrada
switch nargin
    case 2 % La unica operacion admitida es la inversion
        if operacion~='#';error('Dos argumentos y la operacion no es inversion');end
        [fA,cA]=size(A);
        if fA~=cA;error('La matriz a invertir no es cuadrada');end
        if fA~=rank(A);error('La matriz a invertir no es regular');end
        R=inv(A);
    case 3
        [fA,cA]=size(A);[fB,cB]=size(B);
        switch operacion
            case {'+','-'}
                if fA~=fB|cA~=cB;error('Las matrices no son congruentes');end
                if operacion=='-';signo=-1;else;signo=1;end
                R=A+signo*B;
            otherwise
                error('Operación no valida');
        end
    otherwise
        error('Número de argumentos no válido');
end
% ##################### fin del fichero ####################
%
% * * * EJERCICIO 1 * * *
% Añadir las siguientes instrucciones case
            case '*'
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A*B;
            case '\';
                if fA~=cA;error('La matriz a invertir no es cuadrada');end
                if fA~=rank(A);error('La matriz a invertir no es regular');end
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A\B;
            case '/';
                if fB~=cB;error('La matriz a invertir no es cuadrada');end
                if fB~=rank(B);error('La matriz a invertir no es regular');end
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A\B;
% ============
% Apartado 1.2
% ============
%
%
A=[ 112,84,48,28; 84,81,77,-7; 48,77,146,-83; 28,-7,-83,85];
b=[-204;-186;-193;38];
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)
%
% * * * EJERCICIO 2 * * *
%
% Generacion de las matrices
%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
%

SE HA DEJADO INTENCIONADAMENTE EN BLANCO PARA SU REALIZACION
% ##################### 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 * * *
%
disp('Jacobi');K=Kj;
disp(['Normas 1, 2 e Inf:',num2str([norm(K,1),norm(K,2),norm(K,Inf)],'% g']);
disp(['Radio espectral: ',num2str(max(abs(eig(K))),'%g')]);
disp('Gauss-Seidel');K=Kg;
disp(['Normas 1, 2 e Inf:',num2str([norm(K,1),norm(K,2),norm(K,Inf)],'% g']);
disp(['Radio espectral: ',num2str(max(abs(eig(K))),'%g')]);
disp('Relajación');K=Ks;
disp(['Normas 1, 2 e Inf:',num2str([norm(K,1),norm(K,2),norm(K,Inf)],'% g']);
disp(['Radio espectral: ',num2str(max(abs(eig(K))),'%g')]);
% * * * EJERCICIO 4 * * *
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')