|
| From: | Doug Stewart |
| Subject: | Re: bug in residue.m |
| Date: | Fri, 21 Sep 2007 21:30:10 -0500 |
| User-agent: | Thunderbird 1.5.0.13 (Windows/20070809) |
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 address@hidden 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 2The 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 address@hidden <mailto:address@hidden> http://homepage.mac.com/hodelas/tar------------------------------------------------------------------------_______________________________________________ Help-octave mailing list address@hidden https://www.cae.wisc.edu/mailman/listinfo/help-octave_______________________________________________ Help-octave mailing list address@hidden https://www.cae.wisc.edu/mailman/listinfo/help-octave_______________________________________________ Help-octave mailing list address@hidden https://www.cae.wisc.edu/mailman/listinfo/help-octave
| [Prev in Thread] | Current Thread | [Next in Thread] |