function P3

% Very helpful to play with simple small example first with pencil and
% paper. If
%              r   m   m
%       T  =   0   r   m
%              0   0   r
% then it is easy to verify from binomial expansion of  (rI + U)^k that
%
%              r^k     (km)r^(k-1)    (km)r^(k-1) + (k(k-1)m^2)r^(k-2)
%      T^k =    0          r^k                 (km)r^(k-1)
%               0           0                      r^k

% This suggests that convergence to zero will generally take longer if
% (i) r (= rho) is close to one or  (ii) norm(U) (roughly m in the example)
% is large. Regarding the latter, it seems that effect of large U is 
% exaggerated for large n. (For n = 3, norm(U)^2 is involved.)

% Experiments reveal that T^k grows before it is damped out.
% Some random example that look at the effect of rho, norm(U), and n...

close all
figure
set(gcf,'position',[100,100,800,600])
% Vary rho
subplot(3,1,1)
hold on
kMax = 200; kVals = 1:kMax; n = 4; M = 2;
A = MakeT(.90,M,n); vA = PowerValues(A,kMax); 
B = MakeT(.93,M,n); vB = PowerValues(B,kMax); 
C = MakeT(.95,M,n); vC = PowerValues(C,kMax); 
semilogy(kVals,vA,kVals,vB,kVals,vC)
title(sprintf('n = %1d   M = %5.0f',n,M))
ylabel('norm( T^k , 1)')
xlabel('k')
legend('rho = .85','rho = .90','rho = .95')

% Vary M
subplot(3,1,2)
hold on
kMax = 100; kVals = 1:kMax; n = 4; rho = .8;
A = MakeT(rho,200,n); vA = PowerValues(A,kMax); 
B = MakeT(rho,250,n); vB = PowerValues(B,kMax); 
C = MakeT(rho,300,n); vC = PowerValues(C,kMax); 
semilogy(kVals,vA,kVals,vB,kVals,vC)
title(sprintf('n = %1d   rho = %5.2f',n,rho))
ylabel('norm( T^k , 1)')
xlabel('k')
legend('M = 200','M = 250','M = 300')

% Vary n
subplot(3,1,3)
hold on
kMax = 100; kVals = 1:kMax; M = 5; rho = .8;
A = MakeT(rho,M,10);  vA = PowerValues(A,kMax); 
B = MakeT(rho,M,20); vB = PowerValues(B,kMax); 
C = MakeT(rho,M,40); vC = PowerValues(C,kMax); 
semilogy(kVals,vA,kVals,vB,kVals,vC)
title(sprintf('M = %1d   rho = %5.2f',M,rho))
ylabel('norm( T^k , 1)')
xlabel('k')
legend('n = 10','n = 20','n = 40')


function v = PowerValues(T,kMax)
% T is a square matrix and v is a length-kMax row vector with 
% v(k) = norm(T^k,1).
v = zeros(1,kMax);
Tk = T; v(1) = norm(Tk,'fro');
for k=2:kMax
    Tk = T*Tk;
    v(k) = norm(Tk,'fro');
end


function T = MakeT(rho,M,n)
% Generates a random n-by-n upper triangular matrix T with diagonal
% elements = rho and strictly upper triangular part that satisfies
% norm(triu(T,1),'fro') = M. (Any norm will do.)

U = triu(randn(n,n),1); U = M*U/norm(U,'fro');
T = rho*eye(n,n) + U;