help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: residue() confusion


From: Ben Abbott
Subject: Re: residue() confusion
Date: Sun, 23 Sep 2007 07:54:02 -0700 (PDT)

In the case of Matlab the syntax would be

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

For Matlab there is no fourth returned value. Matlab handles mutiplicity by
implication. When multiple roots roots occur, they follow each other
sequentially in the "a" and "p" vectors. Therefore, if the poles are grouped
according to their complex conjugate pairs, the result will not be
compatible with Matlab. See this 
http://www.nabble.com/residue%28%29-confusion-tf4502015.html post  for
Mathwork's description.

So the question is do we want to maintain compatibility with Matlab?


Doug Stewart wrote:
> 
> Thanks for the reply.
> In fact the order that you want is there in the software but I thought 
> it would be nice to have the complex conjugate poles pair up.
> To get Matlab's order just comment out the last 3 lines of code.
> 
> To the List do we want Matlab's exact order?
> Its fine with me ether way.
> 
> Ben Abbott wrote:
>> 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.
>>   
> 
> All the code that I added does not directly touch on the accuracy of the 
> code. I only sorted  the data  a different way before doing the code 
> that  was written years ago.
> 
> By all means if you can make the code better we will all cheer.
> 
> Thanks
> 
> Doug
> 
> 
> 
>> 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
>>>
>>>
>>>     
>>
>>   
> 
> _______________________________________________
> 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#a12847050
Sent from the Octave - General mailing list archive at Nabble.com.



reply via email to

[Prev in Thread] Current Thread [Next in Thread]