|
From: | A. Scottedward Hodel |
Subject: | bug in residue.m |
Date: | Tue, 18 Sep 2007 10:44:24 -0500 |
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 address@hidden |
[Prev in Thread] | Current Thread | [Next in Thread] |