|
»
»
« Return to Thread: bug in residue.m
Re: bug in residue.m
by dastew
::
Rate this Message:
Hi Scott Hodel: I value your thoughts on this fix.
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 back together again
#first the real roots
rp=p(find(imag(p)==0));
rps=sort(rp);
## Now the negative complex roots
ipn=p(find(imag(p)<0));
[xx ii]=sort(abs(imag(ipn)));
ipns=ipn(ii);
## now the positive complex roots
ipp=p(find(imag(p)>0));
[xx ii]=sort(imag(ipp));
ipps=ipp(ii);
## now put all the complex roots together
ips=[ipps' ipns']';
## and put all the roots together
p=[rps' ips']';
## this has sorted the real roots and then the complex roots
## with the complex roots sorted by the imaginary parts in
## such a way that the complex conjugate pair will have the same
## multiplicity for each half of the pair.
and then at the end just before endfunction put this
## now put the complex conjugates back together again.
[p ii]=sort(p);
e=e(ii);
r=r(ii);
This seems to work for all cases that I have tried.
Doug Stewart
Doug Stewart wrote:
> 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
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 back together again
#first the real roots
rp=p(find(imag(p)==0));
rps=sort(rp);
## Now the negative complex roots
ipn=p(find(imag(p)<0));
[xx ii]=sort(abs(imag(ipn)));
ipns=ipn(ii);
## now the positive complex roots
ipp=p(find(imag(p)>0));
[xx ii]=sort(imag(ipp));
ipps=ipp(ii);
## now put all the complex roots together
ips=[ipps' ipns']';
## and put all the roots together
p=[rps' ips']';
## this has sorted the real roots and then the complex roots
## with the complex roots sorted by the imaginary parts in
## such a way that the complex conjugate pair will have the same
## multiplicity for each half of the pair.
and then at the end just before endfunction put this
## now put the complex conjugates back together again.
[p ii]=sort(p);
e=e(ii);
r=r(ii);
This seems to work for all cases that I have tried.
Doug Stewart
Doug Stewart wrote:
> 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
« Return to Thread: bug in residue.m
| Free embeddable forum powered by Nabble | Forum Help |
