# The Matrix Computation Toolbox

### Nick Higham (view profile)

11 Sep 2002 (Updated )

A collection of M-files for carrying out various numerical linear algebra tasks.

ps(A, m, tol, rl, marksize)
```function y = ps(A, m, tol, rl, marksize)
%PS     Dot plot of a pseudospectrum.
%       PS(A, M, TOL, RL) plots an approximation to a pseudospectrum
%       of the square matrix A, using M random perturbations of size TOL.
%       M defaults to a SIZE(A)-dependent value and TOL to 1E-3.
%       RL defines the type of perturbation:
%         RL =  0 (default): absolute complex perturbations of 2-norm TOL.
%         RL =  1:           absolute real perturbations of 2-norm TOL.
%         RL = -1:           componentwise real perturbations of size TOL.
%       The eigenvalues of A are plotted as crosses `x'.
%       PS(A, M, TOL, RL, MARKSIZE) uses the specified marker size instead
%       of a size that depends on the figure size, the matrix order, and M.
%       If MARKSIZE < 0, the plot is suppressed and the plot data is returned
%       as an output argument.
%       PS(A, 0) plots just the eigenvalues of A.

%       For a given TOL, the pseudospectrum of A is the set of
%       pseudo-eigenvalues of A, that is, the set
%       { e : e is an eigenvalue of A+E, for some E with NORM(E) <= TOL }.
%
%       References:
%       L. N. Trefethen, Computation of pseudospectra, Acta Numerica,
%          8:247-295, 1999.
%       L. N. Trefethen, Spectra and pseudospectra, in The Graduate
%          Student's Guide to Numerical Analysis '98, M. Ainsworth,
%          J. Levesley, and M. Marletta, eds., Springer-Verlag, Berlin,
%          1999, pp. 217-250.

if diff(size(A)), error('Matrix must be square.'), end
n = length(A);

if nargin < 5, marksize = 0; end
if nargin < 4, rl = 0; end
if nargin < 3, tol = 1e-3; end
if nargin < 2 | isempty(m), m = 5*max(1, round( 25*exp(-0.047*n) )); end

if m == 0
e = eig(A);
ax = cpltaxes(e);
plot(real(e), imag(e), 'x')
axis(ax), axis('square')
return
end

x = zeros(m*n,1);
i = sqrt(-1);

for j = 1:m
if rl == -1     % Componentwise.
dA = -ones(n) + 2*rand(n);   % Uniform random numbers on [-1,1].
dA = tol * A .* dA;
else
if rl == 0   % Complex absolute.
dA = randn(n) + i*randn(n);
else         % Real absolute.
dA = randn(n);
end
dA = tol/norm(dA)*dA;
end
x((j-1)*n+1:j*n) = eig(A + dA);
end

if marksize >= 0

ax = cpltaxes(x);
h = plot(real(x),imag(x),'.');
axis(ax), axis('square')

% Next block adapted from SPY.M.
if marksize == 0
units = get(gca,'units');
set(gca,'units','points');
pos = get(gca,'position');
nps = 2.4*sqrt(n*m);  % Factor based on number of pseudo-ei'vals plotted.
myguess = round(3*min(pos(3:4))/nps);
marksize = max(1,myguess);
set(gca,'units',units);
end

hold on
e = eig(A);
plot(real(e),imag(e),'x');
set(h,'markersize',marksize);
hold off

else

y = x;

end
```