function amath3132hw2
  disp('Solution for amath 3132 hw #2')
  disp('----------My Information-------------')
  
  Name = 'Jahrul';
  Surname = 'Alam';
  ID = 123456789;
  
  fprintf('  Professor:  %s %s (ID: %d) \n', Name, Surname, ID);
  
  
  disp('  ------------------------------------------------  ')  
  disp('  ')
  
  disp('Hit enter to proceed')
  pause
  clear all; close all;
  
  disp('  ----------   Solution to Question #4  ----------  ')  
  disp('  ')
  disp('4(a)')
  P=[101.3;95.5;89.9;79.5;70.1;61.6;54.0;47.2;35.6;26.4;19.3;14.1;10.3;7.5;5.5;];

  sprintf('L2(p) is %8.2f and verified by matlab %10.5f', matrix_norm(P,2), norm(P,2))
  sprintf('Linf(p) is %8.2f and verified by matlab %10.5f', matrix_norm(P,inf), norm(P,inf))

  rho=[1.225;1.168;1.112;1.007;0.909;0.819;0.736;0.660;0.525;0.413;0.311;0.226;0.165;0.120;0.088;];
  sprintf('L2(rho) is %8.2f and verified by matlab %10.5f', matrix_norm(rho,2), norm(rho,2))
  sprintf('Linf(rho) is %8.2f and verified by matlab %10.5f', matrix_norm(rho,inf), norm(rho,inf))
disp('-----')
  matrix_norm(rho,inf)

  disp('hilbert(n) implements the hilbert matrix')
  disp('4(b)')
  disp('matrix_norm(A,p) computes the matrix norm')
  disp('4(c)');disp('Hit enter to proceed');pause;
  
  % 4(c) L1, L2, L_inf norm of hilbert matrices of size n
  n=25:25:500;
  
  Nrm=ones(length(n),3);
  
  for k=1:length(n)
    A=hilbert(n(k));
    Nrm(k,1)=matrix_norm(A,1);
    Nrm(k,2)=matrix_norm(A,2);
    Nrm(k,3)=matrix_norm(A,inf);
  end
  disp('Norms have been calculated')
  disp('Hit enter to proceed');pause;
  figure(1);clf
  plot(n,Nrm(:,1),'o-',n,Nrm(:,2),'d:',n,Nrm(:,3),'x--')
  set(gca,'fontsize',16)
  legend('L1','L2','L_\infty','Location','East')
  xlabel('n');ylabel('norm')

  disp('Hit enter to proceed');pause;
  disp('5(a)');
  close(1)  
  % 5(a) Jacobi iteration matrix
  A=[4 3 2;0 5 -2;-1 -2 7];
  Jac=Jac_mat(A);
  L1=matrix_norm(Jac,1);
  L2=matrix_norm(Jac,2);
  Linf=matrix_norm(Jac,inf);
  
  disp('----Jac: Jacobi iteration matrix-------')
  disp(Jac);
  fprintf('L1 norm of Jac %d\n',L1);
  fprintf('L2 norm of Jac %d\n',L2);
  fprintf('L_inf norm of Jac %d\n',Linf);
  
  disp('5(b)');disp('Hit enter to proceed');pause;
  % 5(b) solve the system from 3
  A=[1 0 1;-1 1 0;1 2 -3];
  b=[2;0;0];

  x0=[0;0;0];
  fprintf('1st initial guess: x0=\n'); disp(x0)
  [x,iter,res]=jacobi(A,b,x0,1e-5,500);iter
  x0=[10;0;0];
  fprintf('2nd initial guess: x0=\n'); disp(x0)
  [x,iter,res]=jacobi(A,b,x0,1e-5,500);iter
  x0=[0;10;0];
  fprintf('3rd initial guess: x0=\n'); disp(x0)
  [x,iter,res]=jacobi(A,b,x0,1e-5,500);iter
  
  disp('Jacobi iteration method converges slowly...')
  disp('for randomly  selected  reasonable initial guesses')
  Jac=Jac_mat(A);
  L2=max(abs(eig(Jac)));

  fprintf('spect. rad.(Jac) %d\n',L2)  
  fprintf('Large, but <=1 spectral radius => slow convergence\n')
  disp('  ')
  disp('Hit enter to proceed');pause;

  x0=[0;0;0];
  fprintf('1st initial guess: x0=\n'); disp(x0)
  [x,iter,res]=gaussseidel(A,b,x0,1e-5,500);iter
  x0=[10;0;0];
  fprintf('2nd initial guess: x0=\n'); disp(x0)
  [x,iter,res]=gaussseidel(A,b,x0,1e-5,500);iter
  x0=[0;10;0];
  fprintf('3rd initial guess: x0=\n'); disp(x0)
  [x,iter,res]=gaussseidel(A,b,x0,1e-5,500);iter

  disp('Gauss-Seidel iteration method DOES NOT converge...')
  disp('for randomly  selected  reasonable initial guesses')
  Gas=Gas_mat(A);
  L2=max(abs(eig(Gas)));
  
  fprintf('spect. rad.(Gas) %d\n',L2)  
  fprintf('spectral radius=1 => NO convergence\n')
  disp('  ')
  
  disp('6(a)');disp('Hit enter to proceed');pause;
% 6(a)
A=[3 -2 0 0 0;-2 5 1 0 0;0 1 4 -2 0;0 0 -2 5 1;0 0 0 1 2];
b=[5;-9;0;3;0];
x0=0*b;
tols=[1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8];
J_itrs(1:length(tols))=0;
G_itrs(1:length(tols))=0;

%Jacobi method

for k=1:length(tols)
  [x,J_itrs(k),res]=jacobi(A,b,x0,tols(k),200);
end

%Gauss-Seidel method
for k=1:length(tols)
  [x,G_itrs(k),res]=gaussseidel(A,b,x0,tols(k),200);
end
J_itrs
pause
G_itrs
disp('We can use normal plot');
figure(1)
plot(tols,J_itrs,'o-',tols,G_itrs,'d--')
set(gca,'fontsize',16);
xlabel('tolerance');ylabel('#of iterations')
legend('Jacobi','Gauss-Seidel')
disp('Hit enter to see the semi-logorithmic plot')
pause
figure(1);clf
semilogx(tols,J_itrs,'o-',tols,G_itrs,'d--')
set(gca,'fontsize',16);
xlabel('tolerance');ylabel('#of iterations')
legend('Jacobi','Gauss-Seidel')

disp('Hit enter to clean workspace');
pause
close all; clear all;


% 4(a) use: A=hilbert(n)
function A = hilbert(n)
  
  E=ones(n);
  for I=1:n,
    for J=1:n,
      A(I,J)=E(I,J)/(I+J-1);
    end
    end


% 4(b) use: norm=matrix_norm(A,p)
    
function res=matrix_norm(A,p)
  
  if(p==1)
    res=max(sum(A));
  end

  if(p==2)
    sz=size(A);
    if(sz(2)==1)
      res=sqrt(sum(A.^2));
    else
      res=max(abs(eig(A)));
    end
  end

  if(p==inf)
    sz=size(A);
    if(sz(2)==1)
      res=max(A);
    else
      res=max(sum(A'));
    end
  end
  

% 5(a) Jacobi iteration matrix
function res=Jac_mat(A)
  D=diag(diag(A));
  L=tril(A,-1);
  U=triu(A,1);
  
  res=-inv(D)*(L+U);
  
  
% 5(a) G-S iteration matrix
function res=Gas_mat(A)
  D=diag(diag(A));
  L=tril(A,-1);
  U=triu(A,1);
  
  res=-inv(L+D)*(U);
  