% jacobi Script demonstrates Jacobi, Gauss-Seidel and SOR
% 
% The two-point boundary value problem
%    -u'' = 0, 0<x<1
%     u(0) = 1,  u (1) = -1
% with true solution u = 1-2*x is discretized by the difference method
%    -u(j-1) + 2*u(j) - u(j+1) = 0   j = 1:n
% and the discrete equations are solved first by Jacobi's method,
% then by the Gauss-Seidel method with red-black ordering.
% 
n = 20;                        % modify to set the # of mesh points
% 
ux = [1-2*linspace(0,1,n+2)]'; % exact solution
u0 = [1; rand(n,1); -1];       % initial guess: correct BC; random in interior
% 
h   = 1/(n+1);
nd  = 3;                      % power of 10 desired for error
% 
% Jacobi method eigenvalues are {cos (m pi h)}, m=1:n.
% 
u = u0;              % Initial iterate
e   = norm(u-ux);    % Norm of error
rho = cos(pi*h);     % Spectral radius of M_J
% 
% Estimate number of iterations need to reduce
% error to 10^(-nd)
% 
mxit = ceil((nd+log10(e))/(-log10(rho)));
disp(sprintf('JACOBI: Estimate %d iterations required',mxit))
% 
tic                                   % Start a timing clock
for m = 1:mxit
  u(2:n+1) = 0.5*(u(1:n) + u(3:n+2)); % Simultaneous relaxation :-)
  if mod(m,20) == 0                   % Print progress every 20 steps
    e = norm(u-ux);
    disp(sprintf('%4d %10.2e',m,e))
  end
end
disp(sprintf('Time elapsed: %f',toc)) 
% 
% Gauss-Seidel with red-black ordering.
% Eigenvalues are squares of Jacobi eigenvalues.
% 
u = u0;                % Same starting guess
e   = norm(u-ux);      % Error norm
rho = cos(pi*h)^2;     % M_GS spectral radius
% 
% Estimate number of iterations need to reduce
% error to 10^(-nd)
% 
mxit = ceil((nd+log10(e))/(-log10(rho)));
disp(sprintf('GAUSS-SEIDEL: Estimate %d iterations required',mxit))
% 
tic                                       % Start clock
for m = 1:mxit
  u(2:2:n+1) = 0.5*(u(1:2:n)+u(3:2:n+2)); % Update "red" points
  u(3:2:n+1) = 0.5*(u(2:2:n)+u(4:2:n+2)); % Update "black" points
  if mod(m,20) == 0                       % Display progress every 20 steps
    e = norm(u-ux);
    disp(sprintf('%4d %10.2e',m,e))
  end
end
disp(sprintf('Time elapsed: %f',toc))
% 
% SOR, red-black ordering, 
% overrelaxation parameter chosen by user.
% 
u = u0;                % Same starting guess
r = zeros(size(u));    % Allocate space for residual vector
e   = norm(u-ux);      % Error norm
rho = cos(pi*h)^2;     % M_GS spectral radius
% 
w = input('Please enter relaxation parameter 0 < w < 2: ')
% 
% Estimate # of iterations required -
% Solve quadratic equation for SOR eigenvalue corresponding
% to largest Jacobi eigenvalue:
% 
% (lambda + w - 1)^2 = lambda * w^2 * rho_j^2
% 
rhoj = cos(pi*h);
rho = 0.5*(rhoj^2*w^2-2*(w-1)+rhoj*w*sqrt(rhoj^2*w^2-4*(w-1)));
mxit = ceil((nd+log10(e))/(-log10(rho)));
disp(sprintf('SOR(%4.2f): Estimate %d iterations required',w,mxit))
% 
tic                                       % Start clock
for m = 1:mxit
% Update "red" points
  r(2:2:n+1) = u(1:2:n)+u(3:2:n+2)-2*u(2:2:n+1);  % Residual
  u(2:2:n+1) = u(2:2:n+1) + (w/2)*r(2:2:n+1);     % Overrelax 
% Update "black" points
  r(3:2:n+1) = u(2:2:n)+u(4:2:n+2)-2*u(3:2:n+1);  % Residual
  u(3:2:n+1) = u(3:2:n+1) + (w/2)*r(3:2:n+1);     % Overrelax 
  if mod(m,5) == 0                       % Display progress every 5 steps
    e = norm(u-ux);
    disp(sprintf('%4d %10.2e',m,e))
  end
end
disp(sprintf('Time elapsed: %f',toc))