help-octave
[Top][All Lists]
Advanced
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: residue() confusion

From: Doug Stewart
Subject: Re: residue() confusion
Date: Mon, 24 Sep 2007 07:49:53 -0500
User-agent: Thunderbird 1.5.0.13 (Windows/20070809)
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
reply via email to
[Prev in Thread] Current Thread [Next in Thread]