|
»
»
« Return to Thread: bug in residue.m
Re: bug in residue.m
by A. Scottedward Hodel-5
::
Rate this Message:
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
>>
>>
>
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave
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
>>
>>
>
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave
« Return to Thread: bug in residue.m
| Free embeddable forum powered by Nabble | Forum Help |
