Mercurial > hg > octave-lyh
annotate scripts/polynomial/polyeig.m @ 15197:0c0f4eef30f4
* polyeig.m: Missing semicolon in test.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Sat, 18 Aug 2012 18:26:11 -0400 |
parents | ee9b1270c25a |
children | e3dc9ff8e0f2 |
rev | line source |
---|---|
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
1 ## Copyright (C) 2012 Fotios Kasolis |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
2 ## |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
3 ## This file is part of Octave. |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
4 ## |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
5 ## Octave is free software; you can redistribute it and/or modify it |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
6 ## under the terms of the GNU General Public License as published by |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
8 ## your option) any later version. |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
9 ## |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
10 ## Octave is distributed in the hope that it will be useful, but |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
13 ## General Public License for more details. |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
14 ## |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
15 ## You should have received a copy of the GNU General Public License |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
16 ## along with Octave; see the file COPYING. If not, see |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
17 ## <http://www.gnu.org/licenses/>. |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
18 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
19 ## -*- texinfo -*- |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
20 ## @deftypefn {Function File} {@var{z} =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl}) |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
21 ## @deftypefnx {Function File} {[ @var{v}, @var{z} ] =} polyeig (@var{C0}, @var{C1}, @dots{}, @var{Cl}) |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
22 ## |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
23 ## Solve the polynomial eigenvalue problem of degree @var{l}. |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
24 ## |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
25 ## Given a @var{n*n} matrix polynomial @var{C(s)} = @var{C0 + C1 s + |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
26 ## @dots{} + Cl s^l} polyeig solves the eigenvalue problem (@var{C0} + |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
27 ## @var{C1} + @dots{} + @var{Cl})v = 0. Note that the eigenvalues |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
28 ## @var{z} are the zeros of the matrix polynomial. @var{z} is a |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
29 ## @var{lxn} vector and @var{v} is a @var{(n x n)l} matrix with columns |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
30 ## that correspond to the eigenvectors. |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
31 ## |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
32 ## @seealso{eig, eigs, compan} |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
33 ## @end deftypefn |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
34 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
35 ## Author: Fotios Kasolis |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
36 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
37 function [ z, varargout ] = polyeig (varargin) |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
38 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
39 if ( nargout > 2 ) |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
40 print_usage (); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
41 endif |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
42 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
43 nin = numel (varargin); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
44 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
45 n = zeros (1, nin); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
46 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
47 for cnt = 1 : nin |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
48 if ! ( issquare (varargin{cnt}) ) |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
49 error ("polyeig: coefficients must be square matrices"); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
50 endif |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
51 n(cnt) = size (varargin{cnt}, 1); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
52 endfor |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
53 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
54 if numel (unique (n)) > 1 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
55 error ("polyeig: coefficients must have the same dimensions"); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
56 endif |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
57 n = unique (n); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
58 |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
59 ## matrix polynomial degree |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
60 l = nin - 1; |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
61 |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
62 ## form needed matrices |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
63 C = [ zeros(n * (l - 1), n), eye(n * (l - 1)); |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
64 -cell2mat(varargin(1 : end - 1)) ]; |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
65 |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
66 D = [ eye(n * (l - 1)), zeros(n * (l - 1), n); |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
67 zeros(n, n * (l - 1)), varargin{end} ]; |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
68 |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
69 ## solve generalized eigenvalue problem |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
70 if ( isequal (nargout, 1) ) |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
71 z = eig (C, D); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
72 else |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
73 [ z, v ] = eig (C, D); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
74 varargout{1} = v; |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
75 ## return n-element eigenvectors normalized so |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
76 ## that the infinity-norm = 1 |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
77 z = z(1:n,:); |
15190
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
78 ## max() takes the abs if complex: |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
79 t = max (z); |
ee9b1270c25a
polyeig.m: Style fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15189
diff
changeset
|
80 z /= diag (t); |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
81 endif |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
82 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
83 endfunction |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
84 |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
85 %!test |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
86 %! C0 = [8, 0; 0, 4]; C1 = [1, 0; 0, 1]; |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
87 %! [v,z] = polyeig (C0, C1); |
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
88 %! assert (isequal (z(1), -8), true); |
15197
0c0f4eef30f4
* polyeig.m: Missing semicolon in test.
John W. Eaton <jwe@octave.org>
parents:
15190
diff
changeset
|
89 %! d = C0*v + C1*v*z; |
15189
2ad5e6212cd7
polyeig: remove DOS newlines
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15187
diff
changeset
|
90 %! assert (isequal (norm(d), 0.0), true); |