|
2414
|
1 function b=hup(C) |
|
|
2 %HUP(C) tests if the polynomial C is a Hurwitz-Polynomial. |
|
|
3 % It tests if all roots lie in the left half of the complex |
|
|
4 % plane |
|
|
5 % B=hup(C) is the same as |
|
|
6 % B=all(real(roots(c))<0) but much faster. |
|
|
7 % The Algorithm is based on the Routh-Scheme. |
|
|
8 % C are the elements of the Polynomial |
|
|
9 % C(1)*X^N + ... + C(N)*X + C(N+1). |
|
|
10 % |
|
|
11 % HUP2 works also for multiple polynomials, |
|
|
12 % each row a poly - Yet not tested |
|
|
13 % |
|
|
14 % REFERENCES: |
|
|
15 % F. Gausch "Systemtechnik", Textbook, University of Technology Graz, 1993. |
|
|
16 % Ch. Langraf and G. Schneider "Elemente der Regeltechnik", Springer Verlag, 1970. |
|
|
17 |
|
4907
|
18 % $Id$ |
|
12580
|
19 % Copyright (c) 1995-1998,2008 by Alois Schloegl <alois.schloegl@gmail.com> |
|
4907
|
20 % |
|
|
21 % This program is free software: you can redistribute it and/or modify |
|
|
22 % it under the terms of the GNU General Public License as published by |
|
|
23 % the Free Software Foundation, either version 3 of the License, or |
|
|
24 % (at your option) any later version. |
|
2414
|
25 % |
|
4907
|
26 % This program is distributed in the hope that it will be useful, |
|
|
27 % but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
28 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
29 % GNU General Public License for more details. |
|
2414
|
30 % |
|
4907
|
31 % You should have received a copy of the GNU General Public License |
|
|
32 % along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
2414
|
33 |
|
|
34 [lr,lc] = size(c); |
|
|
35 |
|
|
36 % Strip leading zeros and throw away. |
|
|
37 % not considered yet |
|
|
38 %d=(c(:,1)==0); |
|
|
39 |
|
|
40 % Trailing zeros mean there are roots at zero |
|
|
41 b=(c(:,lc)~=0); |
|
|
42 lambda=b; |
|
|
43 |
|
|
44 n=zeros(lc); |
|
|
45 if lc>3 |
|
|
46 n(4:2:lc,2:2:lc-2)=1; |
|
|
47 end; |
|
|
48 while lc>1 |
|
|
49 lambda(b)=c(b,1)./c(b,2); |
|
|
50 b = b & (lambda>=0) & (lambda<Inf); |
|
|
51 c=c(:,2:lc)-lambda(:,ones(1,lc-1)).*(c*n(1:lc,1:lc-1)); |
|
|
52 lc=lc-1; |
|
|
53 end; |
