%MCTDEMO Demonstration of Matrix Computation Toolbox.
% N. J. Higham.
% The Matrix Computation Toolbox contains test matrices, matrix
% factorizations, visualization functions, direct search optimization
% functions, and other miscellaneous functions.
% The version of the toolbox is
% For this demonstration you will need to view both the command window
% and one figure window.
% This demonstration emphasises graphics and shows only
% some of the features of the toolbox.
pause % Press any key to continue after pauses.
% A list of test matrices available in MATLAB and in the Toolbox (all full,
% square, and of arbitrary dimension) is obtained by typing `matrix':
% The FV command plots the boundary of the field of values of a matrix
% (the set of all Rayleigh quotients) and plots the eigenvalues as
% crosses (`x'). Here are some examples:
% Here is the field of values of the 10-by-10 Grcar matrix:
% Next, we form a random orthogonal matrix and look at its field of values.
% The boundary is the convex hull of the eigenvalues since A is normal.
A = gallery('randsvd',10, 1);
% The PS command plots an approximation to a pseudospectrum of A,
% which is the set of complex numbers that are eigenvalues of some
% perturbed matrix A + E, with the norm of E at most epsilon
% (default: epsilon = 1E-3).
% The eigenvalues of A are plotted as crosses (`x').
% Here are some interesting PS plots.
% First, we use the Kahan matrix, a triangular matrix made up of sines and
% cosines. Here is an approximate pseudospectrum of the 10-by-10 matrix:
% Next, a different way of looking at pseudospectra, via norms of
% the resolvent. (The resolvent of A is INV(z*I-A), where z is a complex
% variable). PSCONT gives a color map with a superimposed contour
% plot. Here we specify a region of the complex plane in
% which the 8-by-8 Kahan matrix is interesting to look at.
pscont(gallery('kahan',8), 0, 20, [0.2 1.2 -0.5 0.5]);
% The triw matrix is upper triangular, made up of 1s and -1s:
% Here is a combined surface and contour plot of the resolvent for N = 11.
% Notice how the repeated eigenvalue 1 `sucks in' the resolvent.
pscont(gallery('triw',11), 2, 15, [-2 2 -2 2]);
% The next PSCONT plot is for the companion matrix of the characteristic
% polynomial of the CHEBSPEC matrix:
A = gallery('chebspec',8); C = compan(poly(A));
% The SHOW command shows the +/- pattern of the elements of a matrix, with
% blanks for zero elements:
pscont(C, 2, 20, [-.1 .1 -.1 .1]);
% The following matrix has a pseudospectrum in the form of a limacon.
n = 25; A = gallery('triw',n,1,2) - eye(n);
sub(A, 6) % Leading principal 6-by-6 submatrix of A.
% We can get a visual representation of a matrix using the SEE
% command, which produces subplots with the following layout:
% | MESH(A) SEMILOGY(SVD(A)) |
% | PS(A) FV(A) |
% where PS is the 1e-3 pseudospectrum and FV is the field of values.
% RSCHUR is an upper quasi-triangular matrix:
% Matlab's MAGIC function produces magic squares:
A = magic(5)
% Using the toolbox routine PNORM we can estimate the matrix p-norm
% for any value of p.
[pnorm(A,1) pnorm(A,1.5) pnorm(A,2) pnorm(A,pi) pnorm(A,inf)]
% As this example suggests, the p-norm of a magic square is
% constant for all p!
% GERSH plots Gershgorin disks. Here are some interesting examples.
% GFPP generates matrices for which Gaussian elimination with partial
% pivoting produces a large growth factor.
% Let's find the growth factor RHO for partial pivoting and complete pivoting
% for a bigger matrix:
A = gfpp(20);
[L, U, P, Q, rho] = gep(A,'p'); % Partial pivoting using Toolbox function GEP.
[L, U, P, Q, rho] = gep(A,'c'); % Complete pivoting using Toolbox function GEP.
% As expected, complete pivoting does not produce large growth here.
% Function MATRIX allows test matrices in the Toolbox and MATLAB to be
% accessed by number. The following piece of code steps through all the
% non-Hermitian matrices of arbitrary dimension, setting A to each
% 10-by-10 matrix in turn. It evaluates the 2-norm condition number and the
% ratio of the largest to smallest eigenvalue (in absolute values).
% c = ; e = ; j = 1;
% for i=1:matrix(0)
% % Double on next line avoids bug in MATLAB 6.5 re. i = 35.
% A = double(matrix(i, 10));
% if ~isequal(A,A') % If not Hermitian...
% c1 = cond(A);
% eg = eig(A);
% e1 = max(abs(eg)) / min(abs(eg));
% % Filter out extremely ill-conditioned matrices.
% if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end
c = ; e = ; j = 1;
% Double on next line avoids bug in MATLAB 6.5 re. i = 35.
A = double(matrix(i, 10));
if ~isequal(A,A') % If not Hermitian...
c1 = cond(A);
eg = eig(A);
e1 = max(abs(eg)) / min(abs(eg));
% Filter out extremely ill-conditioned matrices.
if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end
% The following plots confirm that the condition number can be much
% larger than the extremal eigenvalue ratio.
j = max(size(c));
semilogy(1:j, c, 'x', 1:j, e, 'o'), hold on
semilogy(1:j, c, '-', 1:j, e, '--'), hold off
title('cond: x, eig\_ratio: o')
% Finally, here are three interesting pseudospectra based on pentadiagonal
% Toeplitz matrices:
ps(full(gallery('toeppen',32,0,1/2,0,0,1))); % Propeller
ps(inv(full(gallery('toeppen',32,0,1,1,0,.25)))); % Man in the moon
ps(full(gallery('toeppen',32,0,1/2,1,1,1))); % Fish
clear A L U P Q V D