Re: bug in residue.m

[Henry F. Mollet]

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 address@hidden 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 address@hidden 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 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
>> 
>>   
>

Figure 1.pdf
Description: Binary data