Re: residue() confusion

[Ben Abbott]

I'd like to suggest that we also incorporate the reciprocity that is part of
Matlab's version of residue.m

Specifically,

[r,p,k] = residue(num,den);

[num,den] = residue(r,p,k);

To accomplish this in a consistent manner a common routine is needed to
determine pole multiplicity.

In the case of Matlab, that function is called mpoles.m

I've written one myself, as well as a reciprocal function to return
[num,den] from r,p,k. For the moment, I've called the reciprocal function
rresidue.m

I've attached both. I haven't looked into the detail of either your code or
Doug's to see how easliy the routine mpoles.m might be used to ensure that
both reciprocal parts handle multiplicity in a consistent manner.

Are either you or Doug interested in pursuing this effort?

If the files don't show up on your end, try nabble's forum
http://www.nabble.com/Octave---General-f1897.html

http://www.nabble.com/file/p12891155/rresidue.m rresidue.m 
http://www.nabble.com/file/p12891155/mpoles.m mpoles.m 

A S Hodel-2 wrote:
> 
> I have a kludgy but I think functional fix to residue.m.  It may  
> still get confused if there's a big cluster of poles close by each  
> other,.
> 
> Here's a simple test  code:
> 
>   num = [1 2 3 4]
>   den = conv([1,3*j],1,-3*j])
> den = conv([1,3*j],[1,-3*j])
> den = conv(den,den)
> den = conv(den,[1,2,1])
> [r,p,m,e] = residue(num,den)
> 
> with output:
> 
> r =
> 
>     0.0280000 + 0.0000000i
>     0.0200000 - 0.0000000i
>    -0.0140000 + 0.0017037i
>    -0.0011111 + 0.0633333i
>    -0.0140000 - 0.0017037i
>    -0.0011111 - 0.0633333i
> 
> p =
> 
>    -1.00000 + 0.00000i
>    -1.00000 + 0.00000i
>    -0.00000 - 3.00000i
>    -0.00000 - 3.00000i
>    -0.00000 + 3.00000i
>    -0.00000 + 3.00000i
> 
> m = [](0x0)
> e =
> 
>     1
>     2
>     1
>     2
>     1
>     2
> 
> The patch is as follows (to octave-2.9.14 residue.m)
> 
> diff -c  /usr/local/share/octave/2.9.14/m/polynomial/residue.m .
> *** /usr/local/share/octave/2.9.14/m/polynomial/residue.m       Tue  
> Sep 25 16:30:36 2007
> --- ./residue.m Tue Sep 25 17:24:47 2007
> ***************
> *** 213,218 ****
> --- 213,233 ----
>      index = (abs (p) >= toler & (abs (imag (p)) ./ abs (p) < toler));
>      p(index) = real (p(index));
> 
> +   # sort poles so that multiplicity loop will work
> +
> +   kk = 1;
> +   while(kk < length(p))
> +     cp = p(kk);  % current pole
> +     idx = find( abs(p - cp) < toler );  % find poles close to this one
> +     if(length(idx) > 1 )  % if multiplicity
> +       oidx = find(abs(p - cp) >= toler);  % get the rest of the poles
> +       mp = p(idx);   % get multiple poles
> +       % reorder poles and set these poles equal.
> +       p = [cp*ones(length(idx),1); p(oidx)];
> +       kk += length(idx);
> +     endif
> +   endwhile
> +
>      ## Find the direct term if there is one.
> 
>      if (lb >= la)
> 
> A. Scottedward Hodel address@hidden
> http://homepage.mac.com/hodelas/tar
> 
> 
> On Sep 22, 2007, at 6:17 PM, Henry F. Mollet wrote:
> 
>> Your concern is justified. I don't know how to do partial fractions  
>> by hand
>> when there is multiplicity. Therefore I checked results by hand  
>> using s =
>> linspace (-4i, 4i, 9) as a first check. It appears that Matlab  
>> results are
>> correct if I take into account multiplicity of [1 2 1 2]. Octave  
>> results
>> appear to be incorrect.
>> Henry
>> octave-2.9.14:29> s =
>>   -0 - 4i   0 - 3i   0 - 2i   0 - 1i   0 + 0i   0 + 1i   0 + 2i   0  
>> + 3i   0
>> + 4i
>>
>> Using left hand side of equation:
>> octave-2.9.14:30> y=(s.^2 + 1)./(s.^4 + 18*s.^2 + 81)
>> y =
>>  Columns 1 through 6:
>>   -0.30612 + 0.00000i       NaN +     NaNi  -0.12000 - 0.00000i    
>> 0.00000 -
>> 0.00000i   0.01235 - 0.00000i   0.00000 - 0.00000i
>>  Columns 7 through 9:
>>   -0.12000 + 0.00000i       NaN +     NaNi  -0.30612 + 0.00000i
>>
>> Using right hand side of equation with partial fraction given by  
>> Matlab:
>> octave-2.9.14:31> yMatlab= (0 - 0.0926i)./(s-3i) + (0.2222 -
>> 0.0000i)./(s-3i).^2 + (0 + 0.0926i)./(s+3i) + (0.2222 + 0.0000i)./(s 
>> +3i).^2
>>
>> yMatlab =
>>  Columns 1 through 6:
>>   -0.30611 + 0.00000i       NaN +     NaNi  -0.11997 + 0.00000i    
>> 0.00001 +
>> 0.00000i   0.01236 + 0.00000i   0.00001 + 0.00000i
>>  Columns 7 through 9:
>>   -0.11997 + 0.00000i       NaN +     NaNi  -0.30611 + 0.00000i
>>
>> Using right hand side of equation with partial fraction given by  
>> Octave:
>> octave-2.9.14:32> yOctave=(-3.0108e+06 - 1.9734e+06i)./(s-3i) +  
>> (3.0108e+06
>> + 1.9734e+06i)./(s-3i).^2 + (-3.0108e+06 + 1.9734e+06i)./(3+3i) +
>> (3.0108e+06 - 1.9734e+06i)./(s+3i).^2
>>
>> yOctave =
>>  Columns 1 through 5:
>>   -2.9632e+06 + 2.3337e+06i         NaN +       NaNi  -2.9095e+06 +
>> 2.1230e+06i  -6.2042e+05 + 4.4801e+05i  -1.8417e+05 - 1.7290e+05i
>>  Columns 6 through 9:
>>   -1.2708e+05 - 1.0447e+06i  -1.3307e+06 - 4.0746e+06i         NaN +
>> NaNi  -5.2185e+06 + 1.9084e+06i
>>
>> **********************************
>>
>> on 9/22/07 2:14 PM, Ben Abbott at address@hidden wrote:
>>
>>> I was more concerned about the differences in "a"
>>>
>>> I suppose I'll need to do a derivation and check the correct answer.
>>>
>>> On Sep 22, 2007, at 5:05 PM, Henry F. Mollet wrote:
>>>
>>>> The result for e should be [1 2 1 2] (multiplicity for both poles).
>>>> Note
>>>> that Matlab does not even give e.  My mis-understanding of the
>>>> problem was
>>>> pointed out by Doug Stewart. Doug posted new code yesterday, which
>>>> I've
>>>> tried unsuccessfully, but I cannot be sure that I've implemented
>>>> residual.m
>>>> correctly. The corrected code still produced e = [1 1 1 1] for me.
>>>> Henry
>>>>
>>>>
>>>> on 9/22/07 1:31 PM, Ben Abbott at address@hidden wrote:
>>>>
>>>>>
>>>>> As a result of reading through Hodel's
>>>>> http://www.nabble.com/bug-in-residue.m-tf4475396.html post  I
>>>>> decided to
>>>>> check to see how my Octave installation and my Matlab installation
>>>>> responded
>>>>> to the example
>>>>>
>>>>> Using Matlab v7.3
>>>>> --------------------------
>>>>>  num = [1 0 1];
>>>>>  den = [1 0 18 0 81];
>>>>>  [a,p,k] = residue(num,den)
>>>>>
>>>>> a =
>>>>>
>>>>>         0 - 0.0926i
>>>>>    0.2222 - 0.0000i
>>>>>         0 + 0.0926i
>>>>>    0.2222 + 0.0000i
>>>>>
>>>>>
>>>>> p =
>>>>>
>>>>>    0.0000 + 3.0000i
>>>>>    0.0000 + 3.0000i
>>>>>    0.0000 - 3.0000i
>>>>>    0.0000 - 3.0000i
>>>>>
>>>>>
>>>>> k =
>>>>>
>>>>>      []
>>>>> --------------------------
>>>>>
>>>>> Using Octave 2.9.13 (via Fink) on Mac OSX
>>>>> --------------------------
>>>>>  num = [1 0 1];
>>>>>  den = [1 0 18 0 81];
>>>>>  [a,p,k] = residue(num,den)
>>>>>
>>>>> a =
>>>>>
>>>>>   -3.0108e+06 - 1.9734e+06i
>>>>>   -3.0108e+06 + 1.9734e+06i
>>>>>   3.0108e+06 + 1.9734e+06i
>>>>>   3.0108e+06 - 1.9734e+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
>>>>> --------------------------
>>>>>
>>>>> These are different from both the result that
>>>>> http://www.nabble.com/bug-in-residue.m-tf4475396.html Hodel
>>>>> obtained , as
>>>>> well as different from
>>>>> http://www.nabble.com/bug-in-residue.m-tf4475396.html Mollet's
>>>>>
>>>>> Thoughts anyone?
>>>>>
>>>>
>>>>
>>>
>>
>>
>> _______________________________________________
>> 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
> 
> 

-- 
View this message in context: 
http://www.nabble.com/residue%28%29-confusion-tf4502015.html#a12891155
Sent from the Octave - General mailing list archive at Nabble.com.