--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matrixcomp/ps.m	Wed May 06 14:56:53 2015 +0200
@@ -0,0 +1,90 @@
+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