%% Radios espectrales en JOR, GS-R, SOR
clear
clc
close all

omegas = linspace(1,2,1000);
radios_JOR = 0*omegas;
radios_GS_R = 0*omegas;
radios_SOR = 0*omegas;

A = [3 -1 0 -1; 1 -4 1 1; 0 -1 3 -1; -2 0 -1 4];
b = [-3; 0; 3; 11];


n = size(A,1);

D = diag(diag(A));
E = -tril(A,-1);
F = -triu(A,1);

k = 0;
x = 0*b;
diff = inf;

% Jacobi --> JOR
Q_J =  D \ (E+F);
r_J = D \ b;

% Gauss-Seidel
Q_GS = (D-E) \ F;
r_GS =  (D-E) \ b;

for j = 1:length(omegas)
    w = omegas(j);
    Q_JOR = w*Q_J+(1-w)*eye(n); % JOR
    Q_GS_R = w*Q_GS+(1-w)*eye(n); % GS-R
    Q_SOR = (D-w*E) \ ((1-w)*D + w*F); % SOR
    radios_JOR(j) = abs(eigs(Q_JOR, 1, 'lm'));
    radios_GS_R(j) = abs(eigs(Q_GS_R, 1, 'lm'));
    radios_SOR(j) = abs(eigs(Q_SOR, 1, 'lm'));
end

plot(omegas,radios_JOR, 'k', 'LineWidth',2)
hold on
plot(omegas,radios_GS_R, 'b', 'LineWidth',2)
plot(omegas,radios_SOR, 'r', 'LineWidth',2)
plot(omegas,0*omegas+1, '--k')
xlabel('\omega', 'FontSize',14)
legend('\rho(Q_{JOR})', '\rho(Q_{GS-R})', '\rho(Q_{SOR})', '')

% [~,ind] = min(radios_SOR);
% omegas(ind)
% radios_SOR(ind)

%% probamos con los omegas optimos

clc

x = 0*b;
diff = inf;

w = 1.0861; 
Q_JOR = w*Q_J+(1-w)*eye(n); % JOR (radio optimo = 0.6855)
r_JOR = w*r_J;

w = 1.2873;
Q_GS_R = w*Q_GS+(1-w)*eye(n); % GS-R (radio optimo = 0.3982)
r_GS_R = w*r_GS;

w = 1.1832;
Q_SOR = (D-w*E) \ ((1-w)*D + w*F); % SOR (radio optimo = 0.3312)
r_SOR = w* ((D-w*E) \b);

k = 0;
while diff > 1e-6 && k < 100
    % vec = Q_JOR*x + r_JOR;
    % vec = Q_GS_R*x + r_GS_R;
    vec = Q_SOR*x + r_SOR;
    diff = norm(x-vec)
    x = vec;
    k = k+1;
end

x
k