Assignment 2 Scripts

P1

% P1  Stationary Vectors

clc
home
close all

% Check correctness with a small easy example
disp('Test for correctness...')
disp(' ')
theta = (pi/360)*[0 20 90 150 200 300]
P = SetUp(theta')
[v1,sigma] = StatVec1(theta');
[v2,UCond] = StatVec2(theta');
disp(' ')

disp(sprintf('Norm( ones(1,6)*P - ones(1,6) )  = %15.6e',norm(ones(1,6)*P-ones(1,6))))
disp(sprintf('Smallest Singular Value of I - P = %15.6e',sigma(6)))
disp(sprintf('The condition of U(1:5,1:5)      = %15.6e',UCond))

disp(' ')
disp( ['For StatVec1,   v = ' sprintf('  %9.8f ', v1)])
disp( ['For StatVec2,   v = ' sprintf('  %9.8f ', v2)])


% What happens as two of the points coalesce?

disp(' ')
disp('>>>>>>>>>>>>>>>>>>>>>')
disp(' ')
disp('The v vector as the spacing 1/10^k decreases between two points...')
disp(' ')
disp('       0        1         2         3         4         5         6')
disp('-----------------------------------------------------------------------')
n = 10; 
m = floor(n/2); theta = (pi/180)*linspace(0,360*(n-1)/n,n)';
V = [];
for k=0:6
    phi = theta; phi(m+1) = phi(m) + (theta(m+1)-theta(m))/10^k;
    [v,UCond] = StatVec2(phi);
    V = [V v];  
end
disp(' ')
for i=1:n
    disp(sprintf('  %8.6f',V(i,:)))
end

P2

% P2  Losing Positive Definiteness

% Check correctness with the ill-conditioned 8-by-8 Hilbert matrix
 
  clc
  home
  disp('   k           d(k)         sigma_min(A(1:k,1:k)-d(k)ek*ek'' ')
  disp('------------------------------------------------------------')
  n = 8;
  A = hilb(n);   % A(i,j) = 1/(i+j-1)
  d = PosDefD(A);
  for k = 1:n
     B = A(1:k,1:k);
     B(k,k) = B(k,k) - d(k);
     disp(sprintf('  %2d %20.15f    %20.15f',k,d(k),min(svd(B))))
  end
  
  % Condition Estimation
  
  disp(sprintf('\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n '))
  disp('Examine condEst = norm(A,1)/min(PosDefD(A)) as a 1-norm condition estimator...')
  for n=4:2:10
      A = hilb(n);
      xExact = ones(n,1); b = A*xExact; x = A\b; 
      relErr = norm(x-xExact,1)/norm(xExact,1);
      condEst = norm(A,1)/min(PosDefD(A));
      condExact = cond(A,1);
      disp(' '); disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>'); disp(' ')
      disp(sprintf('A = hilb(%1d), The computed x = ',n)) 
      disp(' '); disp(x)
      disp(sprintf('relative Error = %7.2e\ncondEst*eps    = %7.2e\ncond(A,1)*eps  = %7.2e',...
          relErr,condEst*eps,condExact*eps))
  end

P3

% P3  Vectorizing Multiple Solutions

rand('seed',0); n = 10; m = 100; 
% Generate m well-conditioned n-by-n examples...
A = rand(n,n+1,m);
for p=1:m
    A(1:n,1:n,p) = A(1:n,1:n,p) + n*eye(n,n);
    A(:,:,p) = A(n:-1:1,:,p);
end
% Compare the results 
X = Solve(A);
error = [];
for p=1:m
    error = [error norm(X(:,p) - A(1:n,1:n,p)\A(1:n,n+1,p))];
end
semilogy(error)
title(sprintf('Error = norm( X(:,p) - A(1:n,1:n,p)^{-1} A(1:n,n+1,p) ), p=1:%1d',m))
ylabel('Error')
xlabel('p')



  function x = NoPivSolve(A)
% A is n-by-(n+1) and A(1:n,1:n) is nonsingular
% x solves Ax = A(:,n+1)

  [n,np1] = size(A);
  
% Reduce A to upper triangular and apply updates to A(:,n+1)...
  for k=1:n-1
     for i=k+1:n
        tau = A(i,k)/A(k,k);
        for j=k:n+1
           A(i,j) = A(i,j) - tau*A(k,j);
        end
     end
  end

% Solve the upper triangular system...
  x = zeros(n,1);
  for i = n:-1:1
     x(i) = A(i,n+1);
     for j=i+1:n
        x(i) = x(i) - A(i,j)*x(j);
     end
     x(i) = x(i)/A(i,i);
  end

P4

% P4 Circumscribing Tetrahedrons

clc
home
rand('seed',0);
n = 9;
theta =  -90 + 180*rand(n,1);
phi   = -180 + 360*rand(n,1);
D = Tetra(theta,phi);
disp('The Circumscribing tetrahedrons....')
disp(' ')
disp(' p   q   r   s        Diameter')
disp('----------------------------------')
for p = 1:n
    for q = 1:p-1
        for r = 1:q-1
            for s = 1:r-1
                if D(p,q,r,s) > 0
                    disp(sprintf(' %1d   %1d   %1d   %1d   %15.6e',p,q,r,s,D(p,q,r,s)))
                end
            end
        end
    end
end