%
%   CONV_DEMO_1.M   Study the convergence of the Gauss Seidel method for three
%              different systems (with varying degrees of diagonal dominance)
%
%   This program solves three systems of equations using the Gauss Seidel/SR method 
%   (with the relaxation parameter set at 1.0) as part of a demonstration in the 
%   Lesson 6 Lecture Notes of what we mean by near diagonal dominance.
%
%   This script file calls the SR.M routine to apply the SR method.
%
%   File written by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%

      clear all,  close all  
      format compact
%
%   set parameters for Gauss Seidel method (same for all cases)
      tol = 1e-5;   M = 250;  alf = 1;  xo = zeros(6,1);
%
%   Case 1: fully diagonally dominant
      A = [ 9 3 1 0 0 0; 3 9 3 1 0 0; 1 3 9 3 1 0;
            0 1 3 9 3 1; 0 0 1 3 9 3; 0 0 0 1 3 9];
      b = [10 0 0 0 0 20]';  
      disp('Case 1 solution via x = A\b: ');  x = A\b
      disp('Case 1 solution via Gauss Seidel: ');  [x,k] = sr(A,b,xo,alf,tol,M)
      if k == M, disp(' *** Warning -- Case not converged ***'); end
      L = tril(A,-1);  D = diag(diag(A));  U = triu(A,1);
      B = inv(alf*L + D)*((1-alf)*D - alf*U);   ev = eig(B);
      disp('Case 1 spectral radius: ');   p = max(abs(ev))
%
%   Case 2: "nearly" diagonally dominant
      A = [ 6 3 1 0 0 0; 3 6 3 1 0 0; 1 3 6 3 1 0;
            0 1 3 6 3 1; 0 0 1 3 6 3; 0 0 0 1 3 6];
      b = [10 0 0 0 0 20]';  
      disp('Case 2 solution via x = A\b: ');  x = A\b
      disp('Case 2 solution via Gauss Seidel: ');  [x,k] = sr(A,b,xo,alf,tol,M)
      if k == M, disp(' *** Warning -- Case not converged ***'); end
      L = tril(A,-1);  D = diag(diag(A));  U = triu(A,1);
      B = inv(alf*L + D)*((1-alf)*D - alf*U);   ev = eig(B);
      disp('Case 2 spectral radius: ');   p = max(abs(ev))
%
%   Case 3: "far removed" from being diagonally dominant
      A =  [3 3 1 0 0 0; 3 3 3 1 0 0; 1 3 3 3 1 0;
            0 1 3 3 3 1; 0 0 1 3 3 3; 0 0 0 1 3 3];
      b = [10 0 0 0 0 20]';  
      disp('Case 3 solution via x = A\b: ');  x = A\b
      disp('Case 3 solution via Gauss Seidel: ');  [x,k] = sr(A,b,xo,alf,tol,M)
      if k == M, disp(' *** Warning -- Case not converged ***'); end
      L = tril(A,-1);  D = diag(diag(A));  U = triu(A,1);
      B = inv(alf*L + D)*((1-alf)*D - alf*U);   ev = eig(B);
      disp('Case 3 spectral radius: ');   p = max(abs(ev))
%
%   end of program
 
