Mercurial > octave-nkf

comparison scripts/polynomial/polyeig.m @ 18575:f57148641869 gui-release

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
polyeig.m: Overhaul function for Matlab compatibility (bug #41865). * NEWS: Announce change in return format of polyeig. * polyeig.m: Return a row vector of eigenvalues rather than a matrix with eigenvalues on the diagonal. Improve docstring. Improve input validation. Add %!error tests for input validation.
author Rik <rik@octave.org>
date Sat, 15 Mar 2014 16:29:46 -0700
parents d63878346099
children fe689210525c
comparison
equal deleted inserted replaced
18574:6b4f9cab88d6 18575:f57148641869
24 ## 24 ##
25 ## Given an @var{n*n} matrix polynomial 25 ## Given an @var{n*n} matrix polynomial
26 ## @code{@var{C}(s) = @var{C0} + @var{C1} s + @dots{} + @var{Cl} s^l} 26 ## @code{@var{C}(s) = @var{C0} + @var{C1} s + @dots{} + @var{Cl} s^l}
27 ## polyeig solves the eigenvalue problem 27 ## polyeig solves the eigenvalue problem
28 ## @code{(@var{C0} + @var{C1} + @dots{} + @var{Cl})v = 0}. 28 ## @code{(@var{C0} + @var{C1} + @dots{} + @var{Cl})v = 0}.
29 ##
29 ## Note that the eigenvalues @var{z} are the zeros of the matrix polynomial. 30 ## Note that the eigenvalues @var{z} are the zeros of the matrix polynomial.
30 ## @var{z} is an @var{lxn} vector and @var{v} is an (@var{n} x @var{n})l matrix 31 ## @var{z} is a row vector with @var{n*l} elements. @var{v} is a matrix
31 ## with columns that correspond to the eigenvectors. 32 ## (@var{n} x @var{n}*@var{l}) with columns that correspond to the
33 ## eigenvectors.
32 ## 34 ##
33 ## @seealso{eig, eigs, compan} 35 ## @seealso{eig, eigs, compan}
34 ## @end deftypefn 36 ## @end deftypefn
35 37
36 ## Author: Fotios Kasolis 38 ## Author: Fotios Kasolis
37 39
38 function [ z, varargout ] = polyeig (varargin) 40 function [z, v] = polyeig (varargin)
39 41
40 if ( nargout > 2 ) 42 if (nargin < 1 || nargout > 2)
41 print_usage (); 43 print_usage ();
42 endif 44 endif
43 45
44 nin = numel (varargin); 46 nin = numel (varargin);
47 n = rows (varargin{1});
45 48
46 n = zeros (1, nin); 49 for i = 1 : nin
47 50 if (! issquare (varargin{i}))
48 for cnt = 1 : nin 51 error ("polyeig: coefficients must be square matrices");
49 if (! issquare (varargin{cnt}))
50 error ("polyeig: coefficients must be square matrices");
51 endif 52 endif
52 n(cnt) = size (varargin{cnt}, 1); 53 if (rows (varargin{i}) != n)
54 error ("polyeig: coefficients must have the same dimensions");
55 endif
53 endfor 56 endfor
54
55 if (numel (unique (n)) > 1)
56 error ("polyeig: coefficients must have the same dimensions");
57 endif
58 n = unique (n);
59 57
60 ## matrix polynomial degree 58 ## matrix polynomial degree
61 l = nin - 1; 59 l = nin - 1;
62 60
63 ## form needed matrices 61 ## form needed matrices
64 C = [ zeros(n * (l - 1), n), eye(n * (l - 1)); 62 C = [ zeros(n * (l - 1), n), eye(n * (l - 1));
65 -cell2mat(varargin(1 : end - 1)) ]; 63 -cell2mat(varargin(1:end-1)) ];
66 64
67 D = [ eye(n * (l - 1)), zeros(n * (l - 1), n); 65 D = [ eye(n * (l - 1)), zeros(n * (l - 1), n);
68 zeros(n, n * (l - 1)), varargin{end} ]; 66 zeros(n, n * (l - 1)), varargin{end} ];
69 67
70 ## solve generalized eigenvalue problem 68 ## solve generalized eigenvalue problem
71 if ( isequal (nargout, 1) ) 69 if (nargout == 2)
72 z = eig (C, D); 70 [z, v] = eig (C, D);
73 else 71 v = diag (v);
74 [ z, v ] = eig (C, D); 72 ## return n-element eigenvectors normalized so that the infinity-norm = 1
75 varargout{1} = v;
76 ## return n-element eigenvectors normalized so
77 ## that the infinity-norm = 1
78 z = z(1:n,:); 73 z = z(1:n,:);
79 ## max() takes the abs if complex: 74 ## max() takes the abs if complex:
80 t = max (z); 75 t = max (z);
81 z /= diag (t); 76 z /= diag (t);
77 else
78 z = eig (C, D);
82 endif 79 endif
83 80
84 endfunction 81 endfunction
85 82
86 83
84 %!shared C0, C1
85 %! C0 = [8, 0; 0, 4]; C1 = [1, 0; 0, 1];
86
87 %!test 87 %!test
88 %! C0 = [8, 0; 0, 4]; C1 = [1, 0; 0, 1]; 88 %! z = polyeig (C0, C1);
89 %! assert (z, [-8; -4]);
90
91 %!test
89 %! [v,z] = polyeig (C0, C1); 92 %! [v,z] = polyeig (C0, C1);
90 %! assert (isequal (z(1), -8), true); 93 %! assert (z, [-8; -4]);
94 %! z = diag (z);
91 %! d = C0*v + C1*v*z; 95 %! d = C0*v + C1*v*z;
92 %! assert (isequal (norm(d), 0.0), true); 96 %! assert (norm (d), 0.0);
93 97
98 %% Input validation tests
99 %!error polyeig ()
100 %!error [a,b,c] = polyeig (1)
101 %!error <coefficients must be square matrices> polyeig (ones (3,2))
102 %!error <coefficients must have the same dimensions> polyeig (ones (3,3), ones (2,2))