%------------------------------------------------------------------------- % function [Amin, Amax, Bmin, Bmax] = eigenvalue_plot(A, B); % % Compare the spectrum of two matrices. The eigenvalues are to be % supplied in two arrays A and B, each of dimension n x 2. The first % column contains the real components, the second the imaginary % components. The plots labeled as "half-plots" are constructed by % displaying only the eigenvalues of A with positive imaginary part % and those of B with negative imaginary part. Because the matrices are % assumed to be real, the eigenvalues come in complex conjugate pairs so % it is sufficient to view just the upper or lower half plane for them. % This means only the real-valued eigenvalues will potentially overlay % each other. % % Six plots are created: % % 1. A half-plot comparison on a polar graph, with (x + iy) % converted to (rho, theta). % % 2. A half-plot comparison of the logs of the eigenvalues, % with (x + iy) converted to (rho, theta), and then % (log(rho), theta) is plotted. A log scale is used % for displaying the rho coordinate (distance from the origin). % This allows comparing spectra that have vastly differing % sizes of eigenvalues, and more importantly helps to view % any clustering around the origin. The technique used here, % however, will put the smallest eigenvalue in absolute value % at the origin. Still, this helps show the relative spacings. % % 3. The same as (2), but with (log(rho)-s, theta) graphed where % s is the log of the eigenvalue with smallest norm. Plot 2 % can be deceptive because logarithms of numbers between % zero and one become negative. In that case, 2 shows them % but with the points (log(rho), theta) plotted as % (-log(rho), theta+pi). Shifting by s means the eigenvalues % are actually for a matrix G/s, instead of G. The relative % spacing and relationships are unaffected. % % 4. A standard real vs. imaginary Cartesian plot of all of the % eigenvalues, including conjugate pairs. % % 5. The half-plot version of 4. % % 6. A plot of all of the norms of the eigenvalues. I have not % found a use for this but it may help someone else. % % This file will open a new graphics window to display each plot, % and will optionally return four values: % % Amin : magnitude of smallest eigenvalue of A % Amax : magnitude of largest eigenvalue of A % Bmin : magnitude of smallest eigenvalue of B % Bmax : magnitude of largest eigenvalue of B % % Little effort has been made to skimp on memory or computation. If you % where able to find the complete spectrum of two matrices, then by % comparison the amount of memory and computation eigenvalue_plot() takes % is piddling. % %------------------------------------------------------------------------- % % Randall Bramley % Department of Computer Science % Indiana University % % Started : Tue Aug 22 1995, 10:54:36 EST % Modified : Fri 09 May 2008, 09:53 AM , to include log and norm plots % Last Modified: Mon 09 Mar 2009, 11:35 AM , for A321 % %------------------------------------------------------------------------- function [Amin, Amax, Bmin, Bmax] = eigenvalue_plot(A, B ); %------------------------------------------------------------------------- % % Global parameters that affect the level of output and whether or not to % include legends with the graphs. Generally you'll want: % % GLOBAL_bloviation = logical(0) (false) % GLOBAL_showlegends = logical(1) (true) % % If the user does not set those globals, those are the defaults used. % % Notice the style convention I use is not the one Matlab recommends; % instead of all caps (BLOVIATION) I use the variable with the % string 'GLOBAL_' prefixed. I generally use all capitals for naming % parameters that won't change - e.g., REAL and IMAGINARY below. % Globals tend to be parameters that are routinely changed during the % lifetime of the code. % %------------------------------------------------------------------------- global GLOBAL_showlegends; global GLOBAL_bloviation; if isempty(GLOBAL_showlegends) GLOBAL_showlegends = 1; end if isempty(GLOBAL_bloviation) GLOBAL_bloviation = 1; end %---------------------------------------------------------- % Mnemonic names for real and imaginary parts of A, B %---------------------------------------------------------- REAL = 1; IMAGINARY = 2; [Atheta, Arho] = cart2pol(A(:,REAL), A(:,IMAGINARY)); [Btheta, Brho] = cart2pol(B(:,REAL), B(:,IMAGINARY)); IA = find(Atheta >= 0); IB = find(Btheta <= 0); %----------------------------------------------------------------- % Keep only the ones with positive imaginary component for A, % and keep only the ones with negative imaginary component for B. %----------------------------------------------------------------- Atheta = Atheta(IA); Btheta = Btheta(IB); Arho = Arho(IA); Brho = Brho(IB); Amin = min(Arho); Amax = max(Arho); Bmin = min(Brho); Bmax = max(Brho); n = size(A, 1); %---------------------------------------- % Convert to complex numbers temporarily %---------------------------------------- Ac = A(:,REAL) + i*A(:,IMAGINARY); Bc = B(:,REAL) + i*B(:,IMAGINARY); normsA = zeros(n,1); normsB = zeros(n,1); for k = 1:n normsA(k) = norm(Ac(k)); normsB(k) = norm(Bc(k)); end maxnormsAB = max([normsA; normsB]); minnormsAB = min([normsA; normsB]); %-------------------- % Not needed anymore: %-------------------- clear Ac Bc %------------------------------------------------------------ % Let S be the maximum norm of all the eigenvalues of A and B % to use for setting the axes of the plot. Define s as the % minimum norm. %------------------------------------------------------------ S = maxnormsAB; s = minnormsAB; figure; polar([0], [S], ' '); xlabel('Magnitude'); ylabel('Angle'); %---------------------------------------------------------------- % Another convention to follow: indent blocks between a "hold on" % and the matching "hold off". This helps avoid the problem of % forgetting to turn it off, and have a subsequent plot overwrite % it but with its data not shown properly because of the axes % limits are frozen. %---------------------------------------------------------------- hold on; handleA = polar(Atheta, Arho, 'b+'); handleB = polar(Btheta, Brho, 'ro'); title('Half-Plot Comparison of Eigenvalues of Two Real Matrices'); if GLOBAL_showlegends legend([handleA handleB], '\lambda(A)', ... '\lambda(B)', 'Location', 'Best'); end hold off; lArho = (log(Arho)); logsA = min((lArho)); logSA = max((lArho)); lBrho = (log(Brho)); logsB = min((lBrho)); logSB = max((lBrho)); logs = min([logsA; logsB]); logS = max([logSA; logSB]); if (GLOBAL_bloviation) % Spew out a lot of info with line breaks. disp('---------- In function eigenvalue_plot ----------------- ') disp(' ') message = 'Max norm among B''s eigenvalues is '; disp(sprintf('%s %g', message, max(normsB))); message = 'Min norm among B''s eigenvalues is '; disp(sprintf('%s %g', message, min(normsB))); disp(' ') message = 'Max norm among A''s eigenvalues is '; disp(sprintf('%s %g', message, max(normsA))); message = 'Min norm among A''s eigenvalues is '; disp(sprintf('%s %g', message, min(normsA))); disp(' ') message = 'Max norm among A and B''s eigenvalues is '; disp(sprintf('%s %g', message, S)); message = 'Min norm among A and B''s eigenvalues is '; disp(sprintf('%s %g', message, s)); disp(' ') message = 'Max log of norms of A''s eigenvalues: '; disp(sprintf('%s %g', message, logSA)); message = 'Min log of norms of A''s eigenvalues: '; disp(sprintf('%s %g', message, logsA)); disp(' ') message = 'Max log of norms of B''s eigenvalues: '; disp(sprintf('%s %g', message, logSB)); message = 'Min log of norms of B''s eigenvalues: '; disp(sprintf('%s %g', message, logsB)); disp(' ') message = 'Max log of norms among A and B''s eigenvalues is '; disp(sprintf('%s %g', message, logS)); message = 'Min log of norms among A and B''s eigenvalues is '; disp(sprintf('%s %g', message, logs)); disp(' ') disp('---------- End output from eigenvalue_plot----------------- ') end % ... of verbose spewing. %--------------------------------------------------------------------- % A logarithmic plot, shifted by the min norm (= min rho) encountered. % Because logs can be negative for x in (0,1), subtract s to make % them have all nonnegative entries. %--------------------------------------------------------------------- figure; polar([0], [logS - logs], ' '); %---------------------------------------------------------------- % A minor Matlab idiom here: Getting a multiline label is done by % using a multiline cell array. Remember this trick, it's useful % and commonly needed. %---------------------------------------------------------------- xlabel({sprintf('Magnitudes Divided by exp(%g) ', logs ), ... '(so that min(log(\lambda)) = 0)'}); ylabel('Angle'); hold on; handleA = polar(Atheta, log(Arho) - logs, 'b+'); handleB = polar(Btheta, log(Brho) - logs, 'ro'); title('Shifted Logarithmic Half-Plot of Eigenvalues of Two Real Matrices'); if GLOBAL_showlegends legend([handleA handleB], '\lambda(A)', ... '\lambda(B)', 'Location', 'Best'); end hold off; %---------------------------------------------------------------------- % Less useful are the unshifted logs. The problem here is that this may % give values of rho that are negative, but they get placed on the plot % as if they were (-rho, theta+pi). %---------------------------------------------------------------------- figure; polar([0 0], [logs logS], ' '); xlabel('Unshifted Magnitudes'); ylabel('Angle'); hold on; handleA = polar(Atheta, log(Arho), 'b+'); handleB = polar(Btheta, log(Brho), 'ro'); title('Unshifted Logarithmic Half-Plot of Eigenvalues of Two Real Matrices'); if GLOBAL_showlegends legend([handleA handleB], '\lambda(A)', ... '\lambda(B)', 'Location', 'Best'); end hold off; %------------------------------------------------------- % Cartesian plot of A and B's eigenvalues %------------------------------------------------------- figure; plot(A(:,REAL), A(:,IMAGINARY), 'b+', ... B(:,REAL), B(:,IMAGINARY), 'ro'); title('Full Spectrum of A and B'); xlabel('Real Part') ylabel('Imaginary Part') if GLOBAL_showlegends legend('\lambda(A)', '\lambda(B)'); end %------------------------------------------------------- % Cartesian log plot of A and B (with angle theta) %------------------------------------------------------- figure; plot(1:length(normsA), sort(log(normsA)), 'b+', ... 1:length(normsB), sort(log(normsB)), 'ro'); title('Full Logarithmic Spectrum of A and B'); xlabel('Eigenvalue Number') ylabel('Log of Norms') if GLOBAL_showlegends legend('log(\lambda(A))', 'log(\lambda(B))'); end %------------------------------------------------- % Cartesian plot of norms of A and B's eigenvalues %------------------------------------------------- figure; ma = length(normsA); mb = length(normsB); plot(1:ma, sort(normsA), 'b+', ... 1:ma, sort(normsB), 'ro'); title('Full Set of Norms of the Eigenvalues of A and B'); if GLOBAL_showlegends legend('log(\mid\mid\lambda(A)\mid\mid)', ... 'log(\mid\mid\lambda(B)\mid\mid)'); end unstack; return