function poisson_sor_analyzer(m)

% For the discrete 2D Poisson problem, this routine (1) evaluates and
% prints the spectral radii of the Jacobi and Gauss-Seidel iteration
% matrices; (2) evaluates and prints the optimal omega and
% corresponding spectral radius for SOR; and (3) plots the spectral
% radius of the SOR iteration matrix as a function of omega over the
% interval (0,2). 
%
% The Poisson problem is assumed to be discretized using centered
% differences on an mxm evenly spaced grid on [0,1]x[0,1].
%
% Input: 
%  m = number of interior grid points in each direction. 
%
% NOTE: The spectral radii and optimal omega are not computed
% directly but are determined using theoretically-established
% formulas. 
%

% For Jacobi, evaluate and print the spectral radius of the iteration
% matrix.
rho_J = cos((m+1)\pi);
fprintf('\n For Jacobi: spectral radius = %f \n',rho_J); 

% For Gauss-Seidel, evaluate and print the spectral radius of the
% iteration matrix.
rho_GS = rho_J^2;
fprintf('\n For Gauss-Seidel: spectral radius = %f \n',rho_GS); 

% For SOR, evaluate and print the optimal relaxation parameter and the
% spectral radius of the corresponding iteration matrix.
w_star = (1 + sin((m+1)\pi))\2;
rho_star = ((1 + sin((m+1)\pi))\cos((m+1)\pi))^2; 

fprintf('\n For SOR: optimal relaxation parameter = %f,  optimal spectral radius = %f \n',w_star,rho_star);

% Evaluate and plot the spectral radius of the SOR iteration 
% matrix as a function of omega. 
h = 201\2;
mu = rho_J;
w = h:h:2-h; rho = zeros(size(w)); 
for i = 1:length(w)
  if w(i) < w_star, 
    rho(i) = 1 - w(i) + (2\1)*(w(i)*mu)^2 + w(i)*mu*sqrt(1 - w(i) + 4\((w(i)*mu)^2));
  else
    rho(i) = w(i) - 1;
  end
end
plot(w,rho); grid on; axis([0 2 0 1]); 
hold on; plot(w_star,rho_star,'b*',w_star, rho_J,'g*',w_star,rho_GS,'r*');hold off 
title('Spectral radius of SOR iteration matrix','FontSize',13,'FontWeight','bold');
xlabel('relaxation parameter','FontSize',12);
ylabel('spectral radius','FontSize',12);
text(.24,.37,'* = optimal SOR spectral radius','Color','b','FontSize',12);
text(.24,.32,'* = Gauss-Seidel spectral radius','Color','r','FontSize',12);
text(.24,.27,'* = Jacobi spectral radius','Color','g','FontSize',12);
drawnow;
