#Métodos Numéricos - 2023. IMERL, Fing
#Código de ejemplo para el ejercicio 18 del práctico 2
clear all;
clc;
clf;

A=[-3,4,0;1,-2,3;-1,9,3];
b=[-3,4,2]';
[n,n]=size(A);
D = diag(diag(A));
E = -tril(A,-1);
F = -triu(A,1);

#Verificar que la solucion del sistema es x=[1,0,1]'
x_sol=[1,0,1]';
norm(A*x_sol-b)

#Implementar la resolución del sistema mediante J y GS
#Demostrar que son divergentes
maxIter=10;
tol=1e-3;
x0=zeros(n,1);

#J
k=1;
diff=inf;
x=x0;
errorJ=zeros(maxIter,1);
QJ=D\(E+F);
rJ=D\b;
#Calcular los vaps de QJ
lambdasJ=eig(QJ)

while (maxIter>k) && (diff>tol)
  x_=QJ*x+rJ;
  diff=norm(x_-x);
  x=x_;
  k=k+1;
  errorJ(k)=diff;
endwhile

#GS
k=1;
diff=inf;
x=x0;
errorGS=zeros(maxIter,1);
QGS=(D-E)\F;
rGS=(D-E)\b;

#Calcular los vaps de QGS
lambdasGS=eig(QGS)

while (maxIter>k) && (diff>tol)
  x_=QGS*x+rGS;
  diff=norm(x_-x);
  x=x_;
  k=k+1;
  errorGS(k)=diff;
endwhile

figure(1);
plot(1:maxIter,errorJ,'LineWidth',2);
hold on;
set(gca,'fontsize', 26);
plot(1:maxIter,errorGS,'LineWidth',2);
title('Evolución de la diferencia entre iterados sucesivos con J y GS','FontSize',28);
xlabel('k')
legend('J','GS');
grid on;

#Considerar JR, GSR y SOR para diferentes omegas
omegas=linspace(1e-5,1-1e-5,100);
radios_jr=zeros(100,1);
radios_gsr=zeros(100,1);
radios_sor=zeros(100,1);

#Hallo el parámetro de relajación óptimo para JR
for i=1:100
  w=omegas(i);
  QJR=w*(D\(E+F))+(1-w)*eye(n);
  radios_jr(i)=max(abs(eig(QJR)));
endfor

#Hallo el parámetro de relajación óptimo para GSR
for i=1:100
  w=omegas(i);
  QGSR=w.*((D-E)\F)+(1-w).*eye(n);
  radios_gsr(i)=max(abs(eig(QGSR)));
endfor

#Hallo el parámetro de relajación óptimo para SOR
for i=1:100
  w=omegas(i);
  QSOR=(D-w*E)\((1-w)*D+w*F);
  radios_sor(i)=max(abs(eig(QSOR)));
endfor

#Considerar JR, GSR y SOR para diferentes omegas
omegas=linspace(1e-5,1-1e-5,100);
radios_jr=zeros(100,1);
radios_gsr=zeros(100,1);
radios_sor=zeros(100,1);

#Hallo los parámetros de relajación óptimos para JR, GSR y SOR
for i=1:100
  w=omegas(i);

  #JR
  QJR=w*(D\(E+F))+(1-w)*eye(n);
  radios_jr(i)=max(abs(eig(QJR)));

  #GSR
  QGSR=w.*((D-E)\F)+(1-w).*eye(n);
  radios_gsr(i)=max(abs(eig(QGSR)));

  #SOR
  QSOR=(D-w*E)\((1-w)*D + w*F);
  radios_sor(i)=max(abs(eig(QSOR)));

endfor

figure(2)
plot(omegas,radios_jr,'LineWidth',2);
hold on;
set(gca,'fontsize', 26);
plot(omegas,radios_gsr,'LineWidth',2);
plot(omegas,radios_sor,'LineWidth',2);
line ([0 1], [1 1], "linestyle", "--", "color", "r");
title('Radio espectral en función de w para JR, GSR y SOR','FontSize',28);
legend('JR','GSR','SOR');
xlabel('w');
grid on;

#Comparación de cantidad de iteraciones entre los métodos
#(c/u con su omega óptimo)
[x, ix_jr]=min(radios_jr);
[x, ix_gsr]=min(radios_gsr);
[x, ix_sor]=min(radios_sor);
wJR=omegas(ix_jr);
wGSR=omegas(ix_gsr);
wSOR=omegas(ix_sor);

QJR=wJR*(D\(E+F))+(1-wJR)*eye(n);
rJR=wJR.*D\b;
QGSR=wGSR.*((D-E)\F)+(1-wGSR).*eye(n);
rGSR=wGSR.*(D-E)\b;
QSOR=(D-wSOR*E)\((1-wSOR)*D + wSOR*F);
rSOR=wSOR*(D-wSOR.*E)\b;

maxIter=100;
tol=1e-3;
x0=zeros(n,1);

k=1;
diff=inf;
x=x0;
errJR=zeros(maxIter,1);

#Cantidad de iteraciones de JR
while (maxIter>k) && (diff>tol)
  x_=QJR*x+rJR;
  diff=norm(x_-x);
  x=x_;
  k=k+1;
  errJR(k)=diff;
endwhile

kJR=k;

k=1;
diff=inf;
x=x0;
errGSR=zeros(maxIter,1);
errGSR(1)=norm(x0-x_sol);

#Cantidad de iteraciones de GSR
while (maxIter>k) && (diff>tol)
  x_=QGSR*x+rGSR;
  diff=norm(x_-x);
  x=x_;
  k=k+1;
  errGSR(k)=diff;
endwhile

kGSR=k;

k=1;
diff=inf;
x=x0;
errSOR=zeros(maxIter,1);
errSOR(1)=norm(x0-x_sol);


#Cantidad de iteraciones de SOR
while (maxIter>k) && (diff>tol)
  x_=QSOR*x+rSOR;
  diff=norm(x_-x);
  x=x_;
  k=k+1;
  errSOR(k)=diff;
endwhile

kSOR=k;

disp('Cantidad de iteraciones');
disp(['JR:  ' num2str(kJR)]);
disp(['GSR: ' num2str(kGSR)]);
disp(['SOR: ' num2str(kSOR)]);

figure(3);
plot(errJR,'LineWidth',2,'DisplayName',strcat('JR-',' ',num2str(kJR),' iteraciones'));
hold on;
set(gca,'fontsize', 26);
plot(errGSR,'LineWidth',2,'DisplayName',strcat('GSR-',' ',num2str(kGSR),' iteraciones'));
plot(errSOR,'LineWidth',2,'DisplayName',strcat('SOR-',' ',num2str(kSOR),' iteraciones'));
title('Evolución de la diferencia entre iterados sucesivos para JR, GSR y SOR','FontSize',30);
legend('show');
xlabel('k');
grid on;