% This determines the eigenvalues for the normal modes of a damped spring
% system and illustrates their movement as the damping parameter changes.
% See Example 4.1, pp. 142-145, of the Demmel book.  

% Specify the number of nodes in the system. 
% n = 4 is the number in the text example. 
n = 4;

M = eye(n,n); M(1,1)=2; M(n,n)=2;
beta = .4;
B = beta*eye(n,n);
K = eye(n,n) - diag(ones(n-1,1),1);
K = K + K';
K(n,n) = K(n,n) - 1;

A = [ -M\B -M\K; eye(n,n) zeros(n,n) ];

% Compute the eigenvalues and plot them in Figure 1. 

lambda = eig(A)
figure(1); plot(real(lambda),imag(lambda),'or','MarkerFaceColor','r');
grid on; axis([-.3 0 -2 2]);
nreals = size(find(imag(lambda)==0),1);
title(['beta = ' num2str(beta) ',   ' num2str(nreals) ' real eigenvalues']);
drawnow;pause;

% Change the scale of the real axis in Figure 2. 

figure(2); plot(real(lambda),imag(lambda),'or','MarkerFaceColor','r');
nreals = size(find(imag(lambda)==0),1);
title(['beta = ' num2str(beta) ',   ' num2str(nreals) ' real eigenvalues']);
grid on; axis([-5.2 0.2 -2 2]);drawnow;
pause;

% Show the movement of the eigenvalues as the damping parameter increases. 

for beta = .4:.01:5
   
    B = beta*eye(n,n);
    A = [ -M\B -M\K; eye(n,n) zeros(n,n) ];
    lambda = eig(A);
    figure(2); plot(real(lambda),imag(lambda),'or','MarkerFaceColor','r');
    nreals = size(find(imag(lambda)==0),1);
    title(['beta = ' num2str(beta) ',   ' num2str(nreals) ' real eigenvalues']);
    grid on;axis([-5.2 0.2 -2 2]);
    %pause(.1); 
end

    

