Re: bug in residue.m

>> on 9/18/07 8:44 AM, A. Scottedward Hodel at hodelas@... 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 hodelas@...
>>> http://homepage.mac.com/hodelas/tar
>>>
>>>
>>> _______________________________________________
>>> Help-octave mailing list
>>> Help-octave@...
>>> https://www.cae.wisc.edu/mailman/listinfo/help-octave
>>>    
>>
>>
>>
>> _______________________________________________
>> Help-octave mailing list
>> Help-octave@...
>> https://www.cae.wisc.edu/mailman/listinfo/help-octave
>>
>>  
>