Re: bug in residue.m

> I think the problem can be fixed by changing the line ( in my version it
> is line #254)
>
>
>  if (abs (p (new_p_index) - p (old_p_index)) < toler)
>
> to
>
>  if (abs (p (new_p_index) - conj(p (old_p_index))) < toler)
>
>
> Scott  would you try this and verify that it is the way to fix it.
>
> Thanks
>
> Doug Stewart
>
>
>
>
> A. Scottedward Hodel 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 hodelas@... <mailto: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
>
>