Mercurial > matrix-functions

comparison toolbox/pscont.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function [x, y, z, m] = pscont(A, k, npts, ax, levels)
2 %PSCONT Contours and colour pictures of pseudospectra.
3 % PSCONT(A, K, NPTS, AX, LEVELS) plots LOG10(1/NORM(R(z))),
4 % where R(z) = INV(z*I-A) is the resolvent of the square matrix A,
5 % over an NPTS-by-NPTS grid.
6 % NPTS defaults to a SIZE(A)-dependent value.
7 % The limits are AX(1) and AX(2) on the x-axis and
8 % AX(3) and AX(4) on the y-axis.
9 % If AX is omitted, suitable limits are guessed based on the
10 % eigenvalues of A.
11 % The eigenvalues of A are plotted as crosses `x'.
12 % K determines the type of plot:
13 % K = 0 (default) PCOLOR and CONTOUR
14 % K = 1 PCOLOR only
15 % K = 2 SURFC (SURF and CONTOUR)
16 % K = 3 SURF only
17 % K = 4 CONTOUR only
18 % The contours levels are specified by the vector LEVELS, which
19 % defaults to -10:-1 (recall we are plotting log10 of the data).
20 % Thus, by default, the contour lines trace out the boundaries of
21 % the epsilon pseudospectra for epsilon = 1e-10, ..., 1e-1.
22 % [X, Y, Z, NPTS] = PSCONT(A, ...) returns the plot data X, Y, Z
23 % and the value of NPTS used.
24 %
25 % After calling this function you may want to change the
26 % color map (e.g., type COLORMAP HOT - see HELP COLOR) and the
27 % shading (e.g., type SHADING INTERP - see HELP INTERP).
28 % For an explanation of the term `pseudospectra' see PS.M.
29 % When A is real and the grid is symmetric about the x-axis, this
30 % routine exploits symmetry to halve the computational work.
31
32 % Colour pseduospectral pictures of this type are referred to as
33 % `spectral portraits' by Godunov, Kostin, and colleagues.
34 % References:
35 % V. I. Kostin, Linear algebra algorithms with guaranteed accuracy,
36 % Technical Report TR/PA/93/05, CERFACS, Toulouse, France, 1993.
37 % L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and
38 % G.A. Watson, eds, Numerical Analysis 1991, Proceedings of the 14th
39 % Dundee Conference, vol. 260, Pitman Research Notes in Mathematics,
40 % Longman Scientific and Technical, Essex, UK, 1992, pp. 234-266.
41
42 Areal = ~norm(imag(A),1);
43
44 if nargin < 5, levels = -10:-1; end
45 e = eig(A);
46 if nargin < 4
47 ax = cpltaxes(e);
48 if Areal, ax(3) = -ax(4); end % Make sure region symmetric about x-axis.
49 end
50 n = max(size(A));
51 if nargin < 3, npts = round( min(max(5, sqrt(20^2*10^3/n^3) ), 30)); end
52 if nargin < 2, k = 0; end
53
54 nptsx = npts; nptsy = npts;
55 Ysymmetry = (Areal & ax(3) == -ax(4));
56
57 x = seqa(ax(1), ax(2), npts);
58 y = seqa(ax(3), ax(4), npts);
59 if Ysymmetry % Exploit symmetry about x-axis.
60 nptsy = ceil(npts/2);
61 y1 = y;
62 y = y(1:nptsy);
63 end
64
65 [xx, yy] = meshgrid(x,y);
66 z = xx + sqrt(-1)*yy;
67 I = eye(n);
68 Smin = zeros(nptsy, nptsx);
69
70 for j=1:nptsx
71 for i=1:nptsy
72 Smin(i,j) = min( svd( z(i,j)*I-A ) );
73 end
74 end
75
76 z = log10( Smin + eps );
77 if Ysymmetry
78 z = [z; z(nptsy-rem(npts,2):-1:1,:)];
79 y = y1;
80 end
81
82 if k == 0 | k == 1
83 pcolor(x, y, z); hold on
84 elseif k == 2
85 surfc(x, y, z); hold on
86 elseif k == 3
87 surf(x, y, z); hold on
88 end
89
90 if k == 0 | k == 4
91 contour(x, y, z, levels); hold on
92 end
93
94 if k ~= 2 & k ~= 3
95 if k == 0 | k == 1
96 s = 'k'; % Black.
97 else
98 s = 'w'; % White.
99 end
100 plot(real(e),imag(e),['x' s]);
101 end
102
103 axis('square');
104 hold off