Assignment 4 Scripts

P1

% P1

% Generate a small example
randn('seed',0);
n = 10;
A = randn(n,n); 
T = hess(A+A');             % Random symmetric n-by-n
d = sort(diag(schur(T)));   % Eigenvalues = d(1) < d(2) < ... < d(n)
a = diag(T);             
b = diag(T,-1);


close all
figure
xmin = d(1)-1;
xmax = d(n)+1;
N = 8;
mu = linspace(d(1),d(n),N);
axis([xmin,xmax+5,1,N+1]);
axis off
hold on
for k=1:N
    m = Sturm(a,b,mu(k));
    plot([xmin xmax],[k k],d,k*ones(n,1),'o',mu(k),k,'*')
    text(xmax+1,k,sprintf(' m = %2d',m))
end
title('m = the number of eigenvalues <= mu')
hold off


% Explain Output

clc
L = 10;
disp('   n      m = #Eigs <= norm(A,1))')
disp('---------------------------------------------')
for n = [10 20 100 200 250 300 320 340 360 380 400 500]
   a = linspace(-L,L,n)';
   b = ones(n-1,1);
   mu = L + 2;
   m = Sturm(a,b,mu);
   disp(sprintf('  %4d       %4d',n,m))
end
disp(' ')
disp('If lambda is an eigenvalue of A, then |lambda| <= norm(A,1).')
disp('So why doesn''t m equal to n?')
disp('Answer this question and submit with output from this problem.')

P2

% P2 
% Check MinEig

% Generate a random symmetric tridiagonal matrix

randn('seed',0)
n = 200;
a = randn(n,1);
b = randn(n-1,1);
A = diag(a) + diag(b,1) + diag(b,-1);

clc
disp('  k         tol       SturmCalls      mu             |mu - mu_true|/tol')
disp('--------------------------------------------------------------------------')
for k=20:20:200
    for tol = [.001  .000001  .000000001]
       [mu,SturmCalls] = MinEig(a(1:k),b(1:k-1),tol);
       % Using schur, find the smallest eigenvalue in absolute value...
       d = sort(diag(schur(A(1:k,1:k)))); 
       m = sum(d<=0);             % Number of eigenvalues <= 0
       if m == 0;
           mu_true = d(1);
       elseif m==k
           mu_true = d(k);
       else
           if abs(d(m)) <= abs(d(m+1))
               mu_true = d(m);
           else
               mu_true = d(m+1);
           end
       end
       ScaledError = abs(mu - mu_true)/tol;
       disp(sprintf('%3d    %10.1e       %4d  %18.12f   %18.12f',k,tol,SturmCalls,mu,ScaledError))
   end
end
       

P3

% P3  Tests SchurPlus

randn('seed',0)
clc
disp('  n       lambda_max(A)        lambda_max(A+v*v'')    norm((A+v*v'')x - mu*x,2)/norm(A+v*v'',inf)')
disp('--------------------------------------------------------------------------------------------------')
for n=20:20:200
    % Generate an n-by-n example
    A = randn(n,n); 
    A = A'*A; 
    v = randn(n,1);
    [Q,D,mu,x] = SchurPlus(A,v);
    [maxEig,k] = max(abs(diag(D)));
    residual = norm((A+v*v' - mu*eye(n,n))*x)/norm(A+v*v',inf);
    disp(sprintf('%4d  %20.15f    %20.15f              %10.3e',n,D(k,k),mu,residual))
    
end