Re: residue() confusion

[Ben Abbott]

Doug,

Your version appears to work for the cases I've tried.

However, the order of the residues and poles are not consistent with
Matlab's

Below is produced by "help residue" at the Matlab prompt.

 RESIDUE Partial-fraction expansion (residues).
    [R,P,K] = RESIDUE(B,A) finds the residues, poles and direct term of
    a partial fraction expansion of the ratio of two polynomials B(s)/A(s).
    If there are no multiple roots,
       B(s)       R(1)       R(2)             R(n)
       ----  =  -------- + -------- + ... + -------- + K(s)
       A(s)     s - P(1)   s - P(2)         s - P(n)
    Vectors B and A specify the coefficients of the numerator and
    denominator polynomials in descending powers of s.  The residues
    are returned in the column vector R, the pole locations in column
    vector P, and the direct terms in row vector K.  The number of
    poles is n = length(A)-1 = length(R) = length(P). The direct term
    coefficient vector is empty if length(B) < length(A), otherwise
    length(K) = length(B)-length(A)+1.

    If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the
    expansion includes terms of the form
                 R(j)        R(j+1)                R(j+m-1)
               -------- + ------------   + ... + ------------
               s - P(j)   (s - P(j))^2           (s - P(j))^m

    [B,A] = RESIDUE(R,P,K), with 3 input arguments and 2 output arguments,
    converts the partial fraction expansion back to the polynomials with
    coefficients in B and A.

so the correct ordering (using your result) would be 

a =

  -0.000000000000000 - 0.092592592592593i
   0.222222222810316 + 0.000000001505971i
   0.000000000000000 + 0.092592592592593i
   0.222222222810316 - 0.000000001505971i

p =

   0.000000016264485 + 2.999999993648592i
  -0.000000016264485 + 3.000000006351409i
   0.000000016264485 - 2.999999993648592i
  -0.000000016264485 - 3.000000006351409i

k = []
e =

                  1
                  2
                  1
                  2

Beyond that, it also appears that the Matlab solution is numerically more
reliable.

It's late here ... tomorrow I'll look over your code to see if I can spot
anything to improve the numerical accuracy. Its not really my strong suit,
but it wont hurt for me to look.


Doug Stewart wrote:
> 
> Here is my results for this question
> 
> num = [1 0 1];
>  den = [1 0 18 0 81]; 
>  [a,p,k,e] = residue(num,den)
> 
> a =
> 
>   1.0e+06  *
> 
>    5.927582766769742 + 2.314767228467131i
>    5.927582730209724 - 2.314767214190160i
>   -5.927582717669299 - 2.314767374340102i
>   -5.927582681109279 + 2.314767360063129i
> 
> p =
> 
>    0.000000016264485 + 2.999999993648592i
>    0.000000016264485 - 2.999999993648592i
>   -0.000000016264485 + 3.000000006351409i
>   -0.000000016264485 - 3.000000006351409i
> 
> k = []
> e =
> 
>                   1
>                   1
>                   1
>                   1
> 
>  >>   [a,p,k,e] = residue2(num,den)
> a =
> 
>   -0.000000000000000 - 0.092592592592593i
>    0.000000000000000 + 0.092592592592593i
>    0.222222222810316 + 0.000000001505971i
>    0.222222222810316 - 0.000000001505971i
> 
> p =
> 
>    0.000000016264485 + 2.999999993648592i
>    0.000000016264485 - 2.999999993648592i
>   -0.000000016264485 + 3.000000006351409i
>   -0.000000016264485 - 3.000000006351409i
> 
> k = []
> e =
> 
>                   1
>                   1
>                   2
>                   2
> 
> 
> 
> the second one is my "fixed" code and this agrees with Matlab.
> 
> You can get my code at:
> 
> www.dougs.homeip.net/octave/residue2.m
> 
> Doug
> 
> 
> 
> 
> 
> 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#a12842668
Sent from the Octave - General mailing list archive at Nabble.com.