help-octave
[Top][All Lists]
Advanced

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

Re: residue() confusion


From: A. Scottedward Hodel
Subject: Re: residue() confusion
Date: Tue, 25 Sep 2007 17:27:48 -0500

I have a kludgy but I think functional fix to residue.m. It may still get confused if there's a big cluster of poles close by each other,.

Here's a simple test  code:

 num = [1 2 3 4]
 den = conv([1,3*j],1,-3*j])
den = conv([1,3*j],[1,-3*j])
den = conv(den,den)
den = conv(den,[1,2,1])
[r,p,m,e] = residue(num,den)

with output:

r =

   0.0280000 + 0.0000000i
   0.0200000 - 0.0000000i
  -0.0140000 + 0.0017037i
  -0.0011111 + 0.0633333i
  -0.0140000 - 0.0017037i
  -0.0011111 - 0.0633333i

p =

  -1.00000 + 0.00000i
  -1.00000 + 0.00000i
  -0.00000 - 3.00000i
  -0.00000 - 3.00000i
  -0.00000 + 3.00000i
  -0.00000 + 3.00000i

m = [](0x0)
e =

   1
   2
   1
   2
   1
   2

The patch is as follows (to octave-2.9.14 residue.m)

diff -c  /usr/local/share/octave/2.9.14/m/polynomial/residue.m .
*** /usr/local/share/octave/2.9.14/m/polynomial/residue.m Tue Sep 25 16:30:36 2007
--- ./residue.m Tue Sep 25 17:24:47 2007
***************
*** 213,218 ****
--- 213,233 ----
    index = (abs (p) >= toler & (abs (imag (p)) ./ abs (p) < toler));
    p(index) = real (p(index));

+   # sort poles so that multiplicity loop will work
+
+   kk = 1;
+   while(kk < length(p))
+     cp = p(kk);  % current pole
+     idx = find( abs(p - cp) < toler );  % find poles close to this one
+     if(length(idx) > 1 )  % if multiplicity
+       oidx = find(abs(p - cp) >= toler);  % get the rest of the poles
+       mp = p(idx);   % get multiple poles
+       % reorder poles and set these poles equal.
+       p = [cp*ones(length(idx),1); p(oidx)];
+       kk += length(idx);
+     endif
+   endwhile
+
    ## Find the direct term if there is one.

    if (lb >= la)

A. Scottedward Hodel address@hidden
http://homepage.mac.com/hodelas/tar


On Sep 22, 2007, at 6:17 PM, 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



reply via email to

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