Mercurial > matrix-functions

comparison toolbox/ps.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function y = ps(A, m, tol, rl, marksize)
2 %PS Dot plot of a pseudospectrum.
3 % PS(A, M, TOL, RL) plots an approximation to a pseudospectrum
4 % of the square matrix A, using M random perturbations of size TOL.
5 % M defaults to a SIZE(A)-dependent value and TOL to 1E-3.
6 % RL defines the type of perturbation:
7 % RL = 0 (default): absolute complex perturbations of 2-norm TOL.
8 % RL = 1: absolute real perturbations of 2-norm TOL.
9 % RL = -1: componentwise real perturbations of size TOL.
10 % The eigenvalues of A are plotted as crosses `x'.
11 % PS(A, M, TOL, RL, MARKSIZE) uses the specified marker size instead
12 % of a size that depends on the figure size, the matrix order, and M.
13 % If MARKSIZE < 0, the plot is suppressed and the plot data is returned
14 % as an output argument.
15 % PS(A, 0) plots just the eigenvalues of A.
16
17 % For a given TOL, the pseudospectrum of A is the set of
18 % pseudo-eigenvalues of A, that is, the set
19 % { e : e is an eigenvalue of A+E, for some E with NORM(E) <= TOL }.
20 %
21 % Reference:
22 % L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and
23 % G.A. Watson, eds, Numerical Analysis 1991, Proceedings of the 14th
24 % Dundee Conference, vol. 260, Pitman Research Notes in Mathematics,
25 % Longman Scientific and Technical, Essex, UK, 1992, pp. 234-266.
26
27 if diff(size(A)), error('Matrix must be square.'), end
28 n = max(size(A));
29
30 if nargin < 5, marksize = 0; end
31 if nargin < 4, rl = 0; end
32 if nargin < 3, tol = 1e-3; end
33 if nargin < 2, m = max(1, round( 25*exp(-0.047*n) )); end
34
35 if m == 0
36 e = eig(A);
37 ax = cpltaxes(e);
38 plot(real(e), imag(e), 'xg')
39 axis(ax);
40 axis('square');
41 return
42 end
43
44 x = zeros(m*n,1);
45 i = sqrt(-1);
46
47 for j = 1:m
48 if rl == -1 % Componentwise.
49 dA = -ones(n) + 2*rand(n); % Uniform random numbers on [-1,1].
50 dA = tol * A .* dA;
51 else
52 if rl == 0 % Complex absolute.
53 dA = randn(n) + i*randn(n);
54 else % Real absolute.
55 dA = randn(n);
56 end
57 dA = tol/norm(dA)*dA;
58 end
59 e = eig(A + dA);
60 x((j-1)*n+1:j*n) = e;
61 end
62
63 if marksize >= 0
64
65 ax = cpltaxes(x);
66 h = plot(real(x),imag(x),'.');
67 axis(ax);
68 axis('square');
69
70 % Next block adapted from SPY.M.
71 if marksize == 0
72 units = get(gca,'units');
73 set(gca,'units','points');
74 pos = get(gca,'position');
75 nps = 2.4*sqrt(n*m); % Factor based on number of pseudo-ei'vals plotted.
76 myguess = round(3*min(pos(3:4))/nps);
77 % [nps myguess]
78 marksize = max(1,myguess);
79 set(gca,'units',units);
80 end
81 set(h,'markersize',marksize);
82 % set(h,'linemarkersize',marksize);
83
84 hold on
85 e = eig(A);
86 plot(real(e),imag(e),'xw');
87 set(h,'markersize',marksize);
88 hold off
89
90 else
91
92 y = x;
93
94 end