Mercurial > matrix-functions
view toolbox/ps.m @ 0:8f23314345f4 draft
Create local repository for matrix toolboxes. Step #0 done.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Wed, 06 May 2015 14:56:53 +0200 |
parents | |
children |
line wrap: on
line source
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 }. % % Reference: % L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and % G.A. Watson, eds, Numerical Analysis 1991, Proceedings of the 14th % Dundee Conference, vol. 260, Pitman Research Notes in Mathematics, % Longman Scientific and Technical, Essex, UK, 1992, pp. 234-266. if diff(size(A)), error('Matrix must be square.'), end n = max(size(A)); if nargin < 5, marksize = 0; end if nargin < 4, rl = 0; end if nargin < 3, tol = 1e-3; end if nargin < 2, m = max(1, round( 25*exp(-0.047*n) )); end if m == 0 e = eig(A); ax = cpltaxes(e); plot(real(e), imag(e), 'xg') 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 e = eig(A + dA); x((j-1)*n+1:j*n) = e; 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); % [nps myguess] marksize = max(1,myguess); set(gca,'units',units); end set(h,'markersize',marksize); % set(h,'linemarkersize',marksize); hold on e = eig(A); plot(real(e),imag(e),'xw'); set(h,'markersize',marksize); hold off else y = x; end