Mercurial > matrix-functions
comparison toolbox/fv.m @ 2:c124219d7bfa draft
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 [f, e] = fv(B, nk, thmax, noplot) | |
| 2 %FV Field of values (or numerical range). | |
| 3 % FV(A, NK, THMAX) evaluates and plots the field of values of the | |
| 4 % NK largest leading principal submatrices of A, using THMAX | |
| 5 % equally spaced angles in the complex plane. | |
| 6 % The defaults are NK = 1 and THMAX = 16. | |
| 7 % (For a `publication quality' picture, set THMAX higher, say 32.) | |
| 8 % The eigenvalues of A are displayed as `x'. | |
| 9 % Alternative usage: [F, E] = FV(A, NK, THMAX, 1) suppresses the | |
| 10 % plot and returns the field of values plot data in F, with A's | |
| 11 % eigenvalues in E. Note that NORM(F,INF) approximates the | |
| 12 % numerical radius, | |
| 13 % max {abs(z): z is in the field of values of A}. | |
| 14 | |
| 15 % Theory: | |
| 16 % Field of values FV(A) = set of all Rayleigh quotients. FV(A) is a | |
| 17 % convex set containing the eigenvalues of A. When A is normal FV(A) is | |
| 18 % the convex hull of the eigenvalues of A (but not vice versa). | |
| 19 % z = x'Ax/(x'x), z' = x'A'x/(x'x) | |
| 20 % => REAL(z) = x'Hx/(x'x), H = (A+A')/2 | |
| 21 % so MIN(EIG(H)) <= REAL(z) <= MAX(EIG(H)) | |
| 22 % with equality for x = corresponding eigenvectors of H. For these x, | |
| 23 % RQ(A,x) is on the boundary of FV(A). | |
| 24 % | |
| 25 % Based on an original routine by A. Ruhe. | |
| 26 % | |
| 27 % References: | |
| 28 % R.A. Horn and C.R. Johnson, Topics in Matrix Analysis, Cambridge | |
| 29 % University Press, 1991, Section 1.5. | |
| 30 % A.S. Householder, The Theory of Matrices in Numerical Analysis, | |
| 31 % Blaisdell, New York, 1964, Section 3.3. | |
| 32 % C.R. Johnson, Numerical determination of the field of values of a | |
| 33 % general complex matrix, SIAM J. Numer. Anal., 15 (1978), | |
| 34 % pp. 595-602. | |
| 35 | |
| 36 if nargin < 2, nk = 1; end | |
| 37 if nargin < 3, thmax = 16; end | |
| 38 thmax = thmax - 1; % Because code below uses thmax + 1 angles. | |
| 39 | |
| 40 iu = sqrt(-1); | |
| 41 [n, p] = size(B); | |
| 42 if n ~= p, error('Matrix must be square.'), end | |
| 43 f = []; | |
| 44 z = zeros(2*thmax+1,1); | |
| 45 e = eig(B); | |
| 46 | |
| 47 % Filter out cases where B is Hermitian or skew-Hermitian, for efficiency. | |
| 48 if norm(skewpart(B),1) == 0 | |
| 49 | |
| 50 f = [min(e) max(e)]; | |
| 51 | |
| 52 elseif norm(symmpart(B),1) == 0 | |
| 53 | |
| 54 e = imag(e); | |
| 55 f = [min(e) max(e)]; | |
| 56 e = iu*e; f = iu*f; | |
| 57 | |
| 58 else | |
| 59 | |
| 60 for m = 1:nk | |
| 61 | |
| 62 ns = n+1-m; | |
| 63 A = B(1:ns, 1:ns); | |
| 64 | |
| 65 for i = 0:thmax | |
| 66 th = i/thmax*pi; | |
| 67 Ath = exp(iu*th)*A; % Rotate A through angle th. | |
| 68 H = 0.5*(Ath + Ath'); % Hermitian part of rotated A. | |
| 69 [X, D] = eig(H); | |
| 70 [lmbh, k] = sort(real(diag(D))); | |
| 71 z(1+i) = rq(A,X(:,k(1))); % RQ's of A corr. to eigenvalues of H | |
| 72 z(1+i+thmax) = rq(A,X(:,k(ns))); % with smallest/largest real part. | |
| 73 end | |
| 74 | |
| 75 f = [f; z]; | |
| 76 | |
| 77 end | |
| 78 % Next line ensures boundary is `joined up' (needed for orthogonal matrices). | |
| 79 f = [f; f(1,:)]; | |
| 80 | |
| 81 end | |
| 82 if thmax == 0; f = e; end | |
| 83 | |
| 84 if nargin < 4 | |
| 85 | |
| 86 ax = cpltaxes(f); | |
| 87 plot(real(f), imag(f),'k') % Plot the field of values | |
| 88 axis(ax); | |
| 89 axis('square'); | |
| 90 | |
| 91 hold on | |
| 92 plot(real(e), imag(e), 'xb') % Plot the eigenvalues too. | |
| 93 hold off | |
| 94 | |
| 95 end |