|
»
»
« Return to Thread: bug in residue.m
Re: bug in residue.m
by dastew
::
Rate this Message:
Hi Scott Hodel:
This seems to work for me. It should be vectorized yet -probably the
find command would shorten it down.
Put this code at line 230 in residue.m
## added by Doug Stewart Sept 2007
## We now separate the real roots from the complex roots
## and sort the real roots
## and sort the complex roots by their imaginary parts
## and then put them bach togethere again
# lp=length(p);
rp(lp)=0;
ip(lp)=0;
r_index=1;
i_index=1;
for k=1:lp
if(isreal(p(k)))
rp(r_index)=p(k);
r_index++;
else
ip(i_index)=p(k);
i_index++;
endif
endfor
rp=rp(1:r_index-1);
rps=sort(rp);
ip=ip(1:i_index-1);
[xx ii]=sort(imag(ip));
ips=ip(ii);
p=[rps ips]'
## this has sorted the real roots and then the complex roots
## with the complex roots sorted by the imaginary parts
Hope this works for all possible roots :-)
Doug Stewart
Doug Stewart wrote:
> I agree.
> I will try and implement that tonight.
> Doug
>
>
> A. Scottedward Hodel wrote:
>
>> 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
>
>
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave
This seems to work for me. It should be vectorized yet -probably the
find command would shorten it down.
Put this code at line 230 in residue.m
## added by Doug Stewart Sept 2007
## We now separate the real roots from the complex roots
## and sort the real roots
## and sort the complex roots by their imaginary parts
## and then put them bach togethere again
# lp=length(p);
rp(lp)=0;
ip(lp)=0;
r_index=1;
i_index=1;
for k=1:lp
if(isreal(p(k)))
rp(r_index)=p(k);
r_index++;
else
ip(i_index)=p(k);
i_index++;
endif
endfor
rp=rp(1:r_index-1);
rps=sort(rp);
ip=ip(1:i_index-1);
[xx ii]=sort(imag(ip));
ips=ip(ii);
p=[rps ips]'
## this has sorted the real roots and then the complex roots
## with the complex roots sorted by the imaginary parts
Hope this works for all possible roots :-)
Doug Stewart
Doug Stewart wrote:
> I agree.
> I will try and implement that tonight.
> Doug
>
>
> A. Scottedward Hodel wrote:
>
>> 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
>
>
_______________________________________________
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 |
