help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: bug in residue.m


From: Doug Stewart
Subject: Re: bug in residue.m
Date: Tue, 18 Sep 2007 21:42:24 -0500
User-agent: Thunderbird 1.5.0.13 (Windows/20070809)

Henry, you missed the difference  in the multiplicity part see  bellow.

Henry F. Mollet wrote:
I had no problem with your example if I understood the problem correctly.
p appeared without having to sort. I had slight problem with inputing "den"
using copy/past from posting.
Octave 2.9.13 on Mac OS X.
Henry

octave-2.9.13:9> num = [1,0,1];
octave-2.9.13:11> den =\302\240 [1,0,18,0,81]
error: invalid character `?' (ASCII 194) near line 11, column 7
parse error:
  syntax error
den =  [1,0,18,0,81]
         ^
octave-2.9.13:11> den = [1,0,18,0,81];

octave-2.9.13:13> [a,p,k,e] = residue(num,den)
a =

  8.4492e+06 - 3.9658e+06i
  8.4492e+06 + 3.9658e+06i
  -8.4492e+06 + 3.9658e+06i
  -8.4492e+06 - 3.9658e+06i

p =

   0.0000 + 3.0000i
   0.0000 - 3.0000i
  -0.0000 + 3.0000i
  -0.0000 - 3.0000i

k = [](0x0)
e =

   1
   1
   1
   1


you have 1
              1
              1
              1


octave-2.9.13:14> help residue
*****************************************************


on 9/18/07 8:44 AM, A. Scottedward Hodel at address@hidden wrote:

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

It should be like above.


I agree with  Hodel

I am going to look at the code too, but I haven had time yet:-(

Doug Stewart
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
http://homepage.mac.com/hodelas/tar


_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave



_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave




reply via email to

[Prev in Thread] Current Thread [Next in Thread]