Re: bug in residue.m

> That's an improvement, but it leads to the same problem (which only
> occurs in pathological cases): if you have a set of poles that are
> sorted with the same sorting measure (e.g., real part, magnitude, or
> complex part) then there is no guarantee that the next pole (or two)
> in the sequence are what you want.
>
> The only safe way that I can think of is:
> For each pole: (say pole k)
>     - for each pole with the same sorting measure (to within tol); say
> there are N of them
>         - check if pole k is the same as pole k+(N-1) (to within tol)  
> If it is, increase the multiplicity of the pair.
>
> It does require a double loop, but since we assume the poles are
> sorted according to some measure it's only necessary to look ahead the
> next few poles, not through the entire list.
>
> A. Scottedward Hodel hodelas@...
> http://homepage.mac.com/hodelas/tar
>
>
> On Sep 19, 2007, at 8:33 AM, Doug Stewart wrote:
>
>> I now see this is not the complete solution.
>>
>> My thoughts so far are this:
>>
>> -  when the roots are real the conj will do no harm
>> - when the roots are complex  they are sorted into complex conjugate
>> pairs
>>      - what we need to do is see if there is a complex conjugate pair
>>            - and then see if the next  two are a complex conjugate pair
>>                 - it they are then
>>                         see if the two pairs are the same
>>                            - if they are then increase the
>> multiplicity of the second pair.
>>
>> Does this look better?
>>
>> Doug Stewart
>>
>>
>> Doug Stewart wrote:
>>> 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
>>>
>>>
>>
>