|
»
»
bug in residue.m
|
View:
New views
10 Messages
—
Rating Filter:
Alert me
|
 |
bug in residue.m
by A. Scottedward Hodel-5
::
Rate this Message:
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@... _______________________________________________ Help-octave mailing list Help-octave@... https://www.cae.wisc.edu/mailman/listinfo/help-octave |
|
|
Re: bug in residue.m
by Henry F. Mollet
::
Rate this Message:
I had no problem with your example if I understood the problem correctly.
p appeared without having to sort. I had slight problem with inputing "den" using copy/past from posting. Octave 2.9.13 on Mac OS X. Henry octave-2.9.13:9> num = [1,0,1]; octave-2.9.13:11> den =\302\240 [1,0,18,0,81] error: invalid character `?' (ASCII 194) near line 11, column 7 parse error: syntax error >>> den = [1,0,18,0,81] ^ octave-2.9.13:11> den = [1,0,18,0,81]; octave-2.9.13:13> [a,p,k,e] = residue(num,den) a = 8.4492e+06 - 3.9658e+06i 8.4492e+06 + 3.9658e+06i -8.4492e+06 + 3.9658e+06i -8.4492e+06 - 3.9658e+06i p = 0.0000 + 3.0000i 0.0000 - 3.0000i -0.0000 + 3.0000i -0.0000 - 3.0000i k = [](0x0) e = 1 1 1 1 octave-2.9.13:14> help residue ***************************************************** on 9/18/07 8:44 AM, A. Scottedward Hodel at hodelas@... 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@... > 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 |
|
|
Re: bug in residue.m
by dastew
::
Rate this Message:
Henry, you missed the difference in the multiplicity part see bellow.
Henry F. Mollet wrote: > I had no problem with your example if I understood the problem correctly. > p appeared without having to sort. I had slight problem with inputing "den" > using copy/past from posting. > Octave 2.9.13 on Mac OS X. > Henry > > octave-2.9.13:9> num = [1,0,1]; > octave-2.9.13:11> den =\302\240 [1,0,18,0,81] > error: invalid character `?' (ASCII 194) near line 11, column 7 > parse error: > syntax error > >>>> den = [1,0,18,0,81] >>>> > ^ > octave-2.9.13:11> den = [1,0,18,0,81]; > > octave-2.9.13:13> [a,p,k,e] = residue(num,den) > a = > > 8.4492e+06 - 3.9658e+06i > 8.4492e+06 + 3.9658e+06i > -8.4492e+06 + 3.9658e+06i > -8.4492e+06 - 3.9658e+06i > > p = > > 0.0000 + 3.0000i > 0.0000 - 3.0000i > -0.0000 + 3.0000i > -0.0000 - 3.0000i > > k = [](0x0) > e = > > 1 > 1 > 1 > 1 > > you have 1 1 1 1 > octave-2.9.13:14> help residue > ***************************************************** > > > on 9/18/07 8:44 AM, A. Scottedward Hodel at hodelas@... 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 >> >> I agree with Hodel I am going to look at the code too, but I haven had time yet:-( Doug Stewart >> 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@... >> 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 |
|
|
Re: bug in residue.m
by dastew
::
Rate this Message:
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 |
|
|
Re: bug in residue.m
by dastew
::
Rate this Message:
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 |
|
|
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 |
|
|
Re: bug in residue.m
by dastew
::
Rate this Message:
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 |
|
|
Re: bug in residue.m
by Henry F. Mollet
::
Rate this Message:
Thanks. I've managed to plot the function Q(s) that determines the poles. It
shows why there is "multiplicity". Can the plotting of imaginary numbers be done more elegantly? I also know that I'm struggling with 'functions'. Henry octave-2.9.13:1> function [Q] = f(s) > Q=s.^4 +18*s.^2 + 81 ## den = [1,0,18,0,81]; given below > end octave-2.9.13:2> s = linspace (-4i, 4i, 100); octave-2.9.13:3> out=f(s); Q = CUT octave-2.9.13:4> plot (s, out) error: octave_base_value::array_value(): wrong type argument `complex matrix' octave-2.9.13:5> imag_s=imag(s); octave-2.9.13:6> real_out=real(out); octave-2.9.13:7> plot (imag_s, real_out); octave-2.9.13:8> axis ([-4,4]) # because x-axis was from -1 to + 1 only. octave-2.9.13:9> Plot in pdf format is attached. on 9/18/07 7:42 PM, Doug Stewart at dastew@... wrote: > Henry, you missed the difference in the multiplicity part see bellow. > > Henry F. Mollet wrote: CUT >> on 9/18/07 8:44 AM, A. Scottedward Hodel at hodelas@... 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 >>> >>> > It should be like above. > > > I agree with Hodel > > I am going to look at the code too, but I haven had time yet:-( > > Doug Stewart >>> 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@... >>> 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 |
|
|
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 |
|
|
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 |
Figure 1.pdf (19K) | Free embeddable forum powered by Nabble | Forum Help |
