[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: bug in residue.m
From: |
Henry F. Mollet |
Subject: |
Re: bug in residue.m |
Date: |
Tue, 18 Sep 2007 15:58:46 -0700 |
User-agent: |
Microsoft-Entourage/11.1.0.040913 |
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
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
>
> 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