Octave 2.9.13 on Mac OS X:
The m-file below reveals a problem in residue.m, in Octave's
polynomial scripts. I started to debug it, but the
code is fairly intricate. The problem is that the code fails to
detect multiple roots.
Consider the case:
octave:7> num = [1,0,1];
octave:7> den = [1,0,18,0,81];
octave:8> [a,p,k,e] = residue(num,den)
fails to detect the multiple poles at +/- j3 on my machine. The
problem appears to be that residue expects the roots to be returned in
a specific order. The problem in this case is resolved by sorting the
poles by their imaginary parts.
octave:9> %sort poles by imaginary part
octave:9> [a,p,k,e] = residue(num,den)
a =
7.3527e-25 + 9.2593e-02i
2.2222e-01 + 2.3902e-09i
-3.6764e-25 - 9.2593e-02i
2.2222e-01 + 2.3902e-09i
p =
-0.0000 - 3.0000i
0.0000 - 3.0000i
0.0000 + 3.0000i
-0.0000 + 3.0000i
k = [](0x0)
e =
1
2
1
2
The change to residue.m is in the following diff: *Note* This will
fix my problem, but it can still break if two pairs of complex poles
have the same imaginary part, e.g., if you have poles at
1+j, 1-j, -1+j, and -1-j,
if they are sorted in order of imaginary part
-1+j,1+j,-1-j, 1-j,
then the code will still fail to detect the multiplicity. The details
of the code are complicated enough that I can't propose a proper fix
right now, but this patch will fix the problem cited above.
*** /sw/share/octave/2.9.13/m/polynomial/residue.m Fri Sep 7
09:44:44 2007
--- residue.m Tue Sep 18 10:38:20 2007
***************
*** 201,207 ****
## Find the poles.
! p = roots (a);
lp = length (p);
## Determine if the poles are (effectively) zero.
--- 201,207 ----
## Find the poles.
! p = sortcom(roots (a), "im");
lp = length (p);
## Determine if the poles are (effectively) zero.
A. Scottedward Hodel address@hidden <mailto:address@hidden>
http://homepage.mac.com/hodelas/tar
------------------------------------------------------------------------
_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave