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: Sun, 23 Sep 2007 06:35:51 -0500
User-agent: Thunderbird 1.5.0.13 (Windows/20070809)

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






reply via email to

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