%
%  lect26.m:  eigenanalysis for the upwind matrix 
%

%  set parameters

N = 16;
lambda = 0.75;
vect = ones(N,1);

%  upwinding matrix 

mtx = diag((1-lambda)*ones(N,1),0) + diag(lambda*ones(N-1,1),-1);
mtx(1,N) = lambda;

%  eigenanalysis:  evl has e-vals & evc has e-vects as cols

[evc, evl] = eig(mtx);

%  shift plots of e-vects

figure(1);clf
subplot(1,2,1)
hold on
axis([1 N 0.5 N+0.5])
for nn = 1:N
  plot(nn+0.3*real(evc(:,nn)/evc(1,nn)),'o')
  plot(nn+0.3*real(evc(:,nn)/evc(1,nn)))
end
title('\bf eigenvectors for N=16 (real part)')
ylabel('\bf w (shifted vertically)')
xlabel('\bf index')

subplot(1,2,2)
hold on
axis([1 N 0.5 N+0.5])
for nn = 1:N
  plot(nn+0.3*imag(evc(:,nn)/evc(1,nn)),'o')
  plot(nn+0.3*imag(evc(:,nn)/evc(1,nn)))
end
title('\bf eigenvectors for N=16 (imag part)')
ylabel('\bf w (shifted vertically)')
xlabel('\bf index')

%  compare e-vals with VN analysis

figure(2);clf
subplot(2,1,1)
hold on
plot(real(diag(evl)),'x')
plot(imag(diag(evl)),'o')
title('\bf matrix eigenvalues')
ylabel('\bf s (real=x, imag=o)')
xlabel('\bf index')

nn = -N/2:1:N/2-1;
alpha_dx = 2*pi*nn/N;

exp_iom = 1 - lambda*(1 - exp(-i*alpha_dx));

subplot(2,1,2)
hold on
plot(real(exp_iom),'x')
plot(imag(exp_iom),'o')
title('\bf eigenvalues from VN analysis')
ylabel('\bf exp^{i \omega dt} (real=x,imag=o)')
xlabel('\bf index')