% spectra of smoothers for multigrid:
% method 1: Gauss-Seidel
% method 2: Jacobi
% method 3: weighted Jacobi, omega = 4/5
% method 4: SOR, omega = 2/(1 + sin(pi*h))

figure(3);
hold on;
grid;
p = 4;
n = 2^p;
h = 1/n;

xi = pi*h : pi*h : (n - 1)*pi*h;
[kx ky] = meshgrid(xi,xi);


% Jacobi
mu = 0.5*(cos(kx) + cos(ky));
plot3(kx(:),ky(:),abs(mu(:)),'.','color','r','Markersize',20);

% weighted Jacobi
si1 = sin(kx*0.5);
si2 = sin(ky*0.5);
om = 4/5;
mu = 1 - om*(si1.^2 + si2.^2);
plot3(kx(:),ky(:),abs(mu(:)),'.','color','m','Markersize',20);

% Gauss - Seidel
mu = 0.25*(cos(kx) + cos(ky)).^2;
mu = mu(:);
plot3(kx(:),ky(:),abs(mu(:)),'.','color','k','Markersize',20);

% SOR
om = 2/(1 + sin(pi*h));
mu = 0.5*(cos(kx) + cos(ky));
for i = 1 : n - 1
    for j = 1 : n - 1
        a = 1; 
        b = 2*(om - 1) - om^2*mu(i,j)^2;
        c = (om - 1 )^2;
        root = roots([a b c]);
        figure(3);
        plot3(xi(i),xi(j),max(abs(root)),'.','color','b','Markersize',20);
        figure(4)
        hold on;
        plot(real(root),imag(root),'.','color','b','Markersize',20);
    end
end
figure(4);
grid;
