1######################################################################## 2## 3## Copyright (C) 2012-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{z} =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl}) 28## @deftypefnx {} {[@var{v}, @var{z}] =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl}) 29## 30## Solve the polynomial eigenvalue problem of degree @var{l}. 31## 32## Given an @var{n}x@var{n} matrix polynomial 33## 34## @code{@var{C}(@var{s}) = @var{C0} + @var{C1} @var{s} + @dots{} + @var{Cl} 35## @var{s}^@var{l}} 36## 37## @code{polyeig} solves the eigenvalue problem 38## 39## @code{(@var{C0} + @var{C1} @var{z} + @dots{} + @var{Cl} @var{z}^@var{l}) 40## @var{v} = 0}. 41## 42## Note that the eigenvalues @var{z} are the zeros of the matrix polynomial. 43## @var{z} is a row vector with @code{@var{n}*@var{l}} elements. @var{v} is a 44## matrix (@var{n} x @var{n}*@var{l}) with columns that correspond to the 45## eigenvectors. 46## 47## @seealso{eig, eigs, compan} 48## @end deftypefn 49 50function [z, v] = polyeig (varargin) 51 52 if (nargin < 1 || nargout > 2) 53 print_usage (); 54 endif 55 56 n = rows (varargin{1}); 57 58 for i = 1 : nargin 59 if (! issquare (varargin{i})) 60 error ("polyeig: coefficients must be square matrices"); 61 endif 62 if (rows (varargin{i}) != n) 63 error ("polyeig: coefficients must have the same dimensions"); 64 endif 65 endfor 66 67 ## matrix polynomial degree 68 l = nargin - 1; 69 70 ## form needed matrices 71 C = [ zeros(n * (l - 1), n), eye(n * (l - 1)); 72 -cell2mat(varargin(1:end-1)) ]; 73 74 D = [ eye(n * (l - 1)), zeros(n * (l - 1), n); 75 zeros(n, n * (l - 1)), varargin{end} ]; 76 77 ## solve generalized eigenvalue problem 78 if (nargout < 2) 79 z = eig (C, D); 80 else 81 [z, v] = eig (C, D); 82 v = diag (v); 83 ## return n-element eigenvectors normalized so that the infinity-norm = 1 84 z = z(1:n,:); 85 t = max (z); # max() takes the abs if complex. 86 z ./= t; 87 endif 88 89endfunction 90 91 92%!shared C0, C1 93%! C0 = [8, 0; 0, 4]; 94%! C1 = [1, 0; 0, 1]; 95 96%!test 97%! z = polyeig (C0, C1); 98%! assert (z, [-8; -4]); 99 100%!test 101%! [v,z] = polyeig (C0, C1); 102%! assert (z, [-8; -4]); 103%! z = diag (z); 104%! d = C0*v + C1*v*z; 105%! assert (norm (d), 0.0); 106 107## Test input validation 108%!error polyeig () 109%!error [a,b,c] = polyeig (1) 110%!error <coefficients must be square matrices> polyeig (ones (3,2)) 111%!error <coefficients must have the same dimensions> 112%! polyeig (ones (3,3), ones (2,2)) 113