Mercurial > matrix-functions

comparison matrixcomp/mctdemo.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 clc
2 format compact
3 echo on
4 %MCTDEMO Demonstration of Matrix Computation Toolbox.
5 % N. J. Higham.
6
7 % The Matrix Computation Toolbox contains test matrices, matrix
8 % factorizations, visualization functions, direct search optimization
9 % functions, and other miscellaneous functions.
10
11 % The version of the toolbox is
12
13 matrix(-1)
14 echo on
15
16 % For this demonstration you will need to view both the command window
17 % and one figure window.
18 % This demonstration emphasises graphics and shows only
19 % some of the features of the toolbox.
20
21 pause % Press any key to continue after pauses.
22
23 % A list of test matrices available in MATLAB and in the Toolbox (all full,
24 % square, and of arbitrary dimension) is obtained by typing `matrix':
25
26 matrix
27
28 pause
29
30 % The FV command plots the boundary of the field of values of a matrix
31 % (the set of all Rayleigh quotients) and plots the eigenvalues as
32 % crosses (`x'). Here are some examples:
33
34 % Here is the field of values of the 10-by-10 Grcar matrix:
35
36 clf
37 fv(gallery('grcar',10));
38 title('fv(gallery(''grcar'',10))')
39
40 pause
41
42 % Next, we form a random orthogonal matrix and look at its field of values.
43 % The boundary is the convex hull of the eigenvalues since A is normal.
44
45 A = gallery('randsvd',10, 1);
46 fv(A);
47 title('fv(gallery(''randsvd'',10, 1))')
48 pause
49
50 % The PS command plots an approximation to a pseudospectrum of A,
51 % which is the set of complex numbers that are eigenvalues of some
52 % perturbed matrix A + E, with the norm of E at most epsilon
53 % (default: epsilon = 1E-3).
54 % The eigenvalues of A are plotted as crosses (`x').
55 % Here are some interesting PS plots.
56
57 % First, we use the Kahan matrix, a triangular matrix made up of sines and
58 % cosines. Here is an approximate pseudospectrum of the 10-by-10 matrix:
59
60 ps(gallery('kahan',10),25);
61 title('ps(gallery(''kahan'',10),25)')
62 pause
63
64 % Next, a different way of looking at pseudospectra, via norms of
65 % the resolvent. (The resolvent of A is INV(z*I-A), where z is a complex
66 % variable). PSCONT gives a color map with a superimposed contour
67 % plot. Here we specify a region of the complex plane in
68 % which the 8-by-8 Kahan matrix is interesting to look at.
69
70 pscont(gallery('kahan',8), 0, 20, [0.2 1.2 -0.5 0.5]);
71 title('pscont(gallery(''kahan'',8))')
72 pause
73
74 % The triw matrix is upper triangular, made up of 1s and -1s:
75
76 gallery('triw',4)
77
78 % Here is a combined surface and contour plot of the resolvent for N = 11.
79 % Notice how the repeated eigenvalue 1 `sucks in' the resolvent.
80
81 pscont(gallery('triw',11), 2, 15, [-2 2 -2 2]);
82 title('pscont(gallery(''triw'',11))')
83 pause
84
85 % The next PSCONT plot is for the companion matrix of the characteristic
86 % polynomial of the CHEBSPEC matrix:
87
88 A = gallery('chebspec',8); C = compan(poly(A));
89
90 % The SHOW command shows the +/- pattern of the elements of a matrix, with
91 % blanks for zero elements:
92
93 show(C)
94
95 pscont(C, 2, 20, [-.1 .1 -.1 .1]);
96 title('pscont(gallery(''chebspec'',8))')
97 pause
98
99 % The following matrix has a pseudospectrum in the form of a limacon.
100
101 n = 25; A = gallery('triw',n,1,2) - eye(n);
102 sub(A, 6) % Leading principal 6-by-6 submatrix of A.
103 ps(A);
104 pause
105
106 % We can get a visual representation of a matrix using the SEE
107 % command, which produces subplots with the following layout:
108 % /---------------------------------\
109 % | MESH(A) SEMILOGY(SVD(A)) |
110 % | PS(A) FV(A) |
111 % \---------------------------------/
112 % where PS is the 1e-3 pseudospectrum and FV is the field of values.
113 % RSCHUR is an upper quasi-triangular matrix:
114
115 see(rschur(16,18));
116
117 pause
118
119 % Matlab's MAGIC function produces magic squares:
120
121 A = magic(5)
122
123 % Using the toolbox routine PNORM we can estimate the matrix p-norm
124 % for any value of p.
125
126 [pnorm(A,1) pnorm(A,1.5) pnorm(A,2) pnorm(A,pi) pnorm(A,inf)]
127
128 % As this example suggests, the p-norm of a magic square is
129 % constant for all p!
130
131 pause
132
133 % GERSH plots Gershgorin disks. Here are some interesting examples.
134 clf
135 gersh(gallery('lesp',12));
136 title('gersh(gallery(''lesp'',12))')
137 pause
138
139 gersh(gallery('hanowa',10));
140 title('gersh(gallery(''hanowa'',10))')
141 pause
142
143 gersh(gallery('ipjfact',6,1));
144 title('gersh(gallery(''ipjfact'',6,1))')
145 pause
146
147 gersh(gallery('smoke',16,1));
148 title('gersh(gallery(''smoke'',16,1))')
149 pause
150
151 % GFPP generates matrices for which Gaussian elimination with partial
152 % pivoting produces a large growth factor.
153
154 gfpp(6)
155 pause
156
157 % Let's find the growth factor RHO for partial pivoting and complete pivoting
158 % for a bigger matrix:
159
160 A = gfpp(20);
161
162 [L, U, P, Q, rho] = gep(A,'p'); % Partial pivoting using Toolbox function GEP.
163 [rho, 2^19]
164
165 [L, U, P, Q, rho] = gep(A,'c'); % Complete pivoting using Toolbox function GEP.
166 rho
167 % As expected, complete pivoting does not produce large growth here.
168 pause
169
170 % Function MATRIX allows test matrices in the Toolbox and MATLAB to be
171 % accessed by number. The following piece of code steps through all the
172 % non-Hermitian matrices of arbitrary dimension, setting A to each
173 % 10-by-10 matrix in turn. It evaluates the 2-norm condition number and the
174 % ratio of the largest to smallest eigenvalue (in absolute values).
175
176 % c = []; e = []; j = 1;
177 % for i=1:matrix(0)
178 % % Double on next line avoids bug in MATLAB 6.5 re. i = 35.
179 % A = double(matrix(i, 10));
180 % if ~isequal(A,A') % If not Hermitian...
181 % c1 = cond(A);
182 % eg = eig(A);
183 % e1 = max(abs(eg)) / min(abs(eg));
184 % % Filter out extremely ill-conditioned matrices.
185 % if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end
186 % end
187 % end
188 echo off
189
190 c = []; e = []; j = 1;
191 for i=1:matrix(0)
192 % Double on next line avoids bug in MATLAB 6.5 re. i = 35.
193 A = double(matrix(i, 10));
194 if ~isequal(A,A') % If not Hermitian...
195 c1 = cond(A);
196 eg = eig(A);
197 e1 = max(abs(eg)) / min(abs(eg));
198 % Filter out extremely ill-conditioned matrices.
199 if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end
200 end
201 end
202 echo on
203
204 % The following plots confirm that the condition number can be much
205 % larger than the extremal eigenvalue ratio.
206 echo off
207 j = max(size(c));
208 subplot(2,1,1)
209 semilogy(1:j, c, 'x', 1:j, e, 'o'), hold on
210 semilogy(1:j, c, '-', 1:j, e, '--'), hold off
211 title('cond: x, eig\_ratio: o')
212 subplot(2,1,2)
213 semilogy(1:j, c./e)
214 title('cond/eig\_ratio')
215 echo on
216 pause
217
218 % Finally, here are three interesting pseudospectra based on pentadiagonal
219 % Toeplitz matrices:
220
221 clf
222 ps(full(gallery('toeppen',32,0,1/2,0,0,1))); % Propeller
223 pause
224
225 ps(inv(full(gallery('toeppen',32,0,1,1,0,.25)))); % Man in the moon
226 pause
227
228 ps(full(gallery('toeppen',32,0,1/2,1,1,1))); % Fish
229 pause
230
231 echo off
232 clear A L U P Q V D
233 format