Re: residue() confusion

[Ben Abbott]

Doug,

I've attached a diff for the version of residue in 2.9.13 ... don't yet have
2.9.14

I have some additional thoughts but will start a fresh thread.
http://www.nabble.com/file/p12870552/residue4.diff residue4.diff 


Doug Stewart wrote:
> 
> The extra e (4th return value) can be ignored and all will work the same 
> as Matlab.
> I agree with you we should return the same order as matlab.
> So here is my currently final version.
> http://dougs.homeip.net/octave/residue4.m
> 
> would someone who has the know how please do a diff against the current 
> residue
> and make sure that the only difference is:
> 
> 
> ## added by Doug Stewart Sept 2007
> ## We now separate the real roots from the complex roots 
> ## and sort the real roots
> ## and sort the complex roots by their imaginary parts
> ## and then put them back together again
> 
> #first the real roots
> rp=p(find(imag(p)==0));
> rps=sort(rp);
> 
> ## Now the negative complex roots
> ipn=p(find(imag(p)<0));
> [xx ii]=sort(abs(imag(ipn)));
>  ipns=ipn(ii);
> 
> ## now the positive complex roots
> ipp=p(find(imag(p)>0));
> [xx ii]=sort(imag(ipp));
>  ipps=ipp(ii);
>  
>  ## now put all the complex roots together
>  ips=[ipps' ipns']';
>  
>  ## and put all the roots together
> p=[rps' ips']';
> ## this has sorted the real roots and then the complex roots
> ## with the complex roots sorted by the imaginary parts in
> ## such a way that the complex conjugate pair will have the same
> ## multiplicity for each half of the pair. 
> 
> 
> 
> then I think this should be put into source forge.
> 
> Can someone do this for me(us).
> Doug Stewart
> 
> 
> 
> 
> Ben Abbott wrote:
>> 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
>>>
>>>
>>>     
>>
>>   
> 
> _______________________________________________
> 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#a12870552
Sent from the Octave - General mailing list archive at Nabble.com.