Re: modified residues() for matlab compatibility

[Ben Abbott]

On Oct 9, 2007, at 2:49 PM, John W. Eaton wrote:

> On  8-Oct-2007, Ben Abbott wrote:
>
> | I did not get the same error. However, I did get an error from the
> | final test script placed on the bottom of the file. Since I was
> | unable to connect to the CVS over the weekend and I did now know if
> | there is a protocol for test-scripts, I'd placed my testing between
> | the doc-strip and remainder of the function.
> |
> | Since having access to the CVS, I've moved the scripts to the bottom
> | and reconciled the duplicate.
> |
> | However, since this is my 1st time using the "test"  
> functionality, it
> | might be wise to double check my work.
> |
> | Meanwhile, there still remains the error you've brought to my
> | attention. The problem appears to be with the order in which the
> | poles/roots are sorted. The script residue.m calls mpoles.m, which
> | uses the internal sort() routine, to sort (reorder) the roots.
> |
> | It appears that your octave reorders the roots differently than  
> mine.
> | I've verified that the mpoles.m checked into the CVS is the same as
> | the one I'm using (however, there was a prior version I had posted
> | that would produce this problem).
> |
> | I haven't yet had the time to mirror the entire CVS to my system. So
> | I am unable to isolate the problem, and be certain of my guess.
> |
> | However, I think we can verify the hypothesis quickly. Using my FInk
> | install of  Octave, I get the following
> |
> |       b = [1, 0, 1];
> |       a = [1, 0, 18, 0, 81];
> |       [r, p, k, e] = residue(b, a)
> |
> | r =
> |
> |     0.00000 - 0.09259i
> |     0.22222 + 0.00000i
> |     0.00000 + 0.09259i
> |     0.22222 - 0.00000i
> |
> | p =
> |
> |     0.0000 + 3.0000i
> |     0.0000 + 3.0000i
> |    -0.0000 - 3.0000i
> |    -0.0000 - 3.0000i
> |
> | k = [](0x0)
> | e =
> |
> |     1
> |     2
> |     1
> |     2
> |
> | On your system, I suspect the poles are ordered as
> |
> | p =
> |
> |     0.0000 - 3.0000i
> |     0.0000 - 3.0000i
> |    -0.0000 + 3.0000i
> |    -0.0000 + 3.0000i
> |
> | i.e. the -3i pair comes ahead of the +3i pair.
> |
> | Can you verify?
>
> Here is what I see with the current sources:
>
>   octave:1> b = [1, 0, 1];
>   octave:2> a = [1, 0, 18, 0, 81];
>   octave:3> [r, p, k, e] = residue(b, a)
>   r =
>
>     -0.00000 + 0.09259i
>      0.22222 - 0.00000i
>     -0.00000 - 0.09259i
>      0.22222 + 0.00000i
>
>   p =
>
>     -0.0000 - 3.0000i
>     -0.0000 - 3.0000i
>      0.0000 + 3.0000i
>      0.0000 + 3.0000i
>
>   k = [](0x0)
>   e =
>
>      1
>      2
>      1
>      2
>
> | If the problem is that simple, then I'll venture a guess that the
> | sort function has been modified since the 2.9.14 release
>
> I don't see any changes that should affect this result and I get the
> sam result if I use the new residue with Octave 2.9.14.
>
> | ... or ...
> | your CPU is giving slightly different answers than mine. In either
> | event, it makes since to include both results in the test-script
> | (since they are each actually correct).
>
> | Assuming my assumptions are correct, I've modified the script. The
> | patch is below. Please let me know if the syntax for creating the
> | patch is not what you prefer.
>
> This patch did not apply for me.  Please use
>
>   cvs diff -u
>
> in the scripts/polynomial directory where you have Octave checked out.
>
> jwe

*sigh* ... sorry John,

in my rush to get the patch to you, I diff'd the wrong file :-(

... I'm not yet able to get the my local copy of the cvs to compile  
due to some dependency issues with my Fink installation. I hope to  
have that problem licked in the next few days.

In any event, the patch below should work. I'm keeping my fingers  
crossed, since I still do not understand why the order of the roots  
is different for our installations. I'll try starting a new thread to  
isolate the problem.

Index: scripts/polynomial/residue.m
===================================================================
RCS file: /cvs/octave/scripts/polynomial/residue.m,v
retrieving revision 1.27
diff -u -r1.27 residue.m
--- scripts/polynomial/residue.m        8 Oct 2007 19:41:28  
-0000       1.27
+++ scripts/polynomial/residue.m        9 Oct 2007 23:25:01 -0000
@@ -122,31 +122,6 @@
 ## Created: June 1994
 ## Adapted-By: jwe

-%!test
-%! b = [1, 1, 1];
-%! a = [1, -5, 8, -4];
-%! [r, p, k, e] = residue (b, a);
-%! assert ((abs (r - [-2; 7; 3]) < 1e-5
-%! && abs (p - [2; 2; 1]) < 1e-7
-%! && isempty (k)
-%! && e == [1; 2; 1]));
-%! k = [1 0];
-%! [b, a] = residue (r, p, k);
-%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
-%! && abs (a - [1, -5, 8, -4]) < 1e-12));
-
-%!test
-%! b = [1, 0, 1];
-%! a = [1, 0, 18, 0, 81];
-%! [r, p, k, e] = residue(b, a);
-%! assert ((abs (54*r - [-5i; 12; 5i; 12]) < 1e-6
-%! && abs (p - [3i; 3i; -3i; -3i]) < 1e-7
-%! && isempty (k)
-%! && e == [1; 2; 1; 2]));
-%! [br, ar] = residue (r, p, k);
-%! assert ((abs (br - b) < 1e-12
-%! && abs (ar - a) < 1e-12));
-
 function [r, p, k, e] = residue (b, a, varargin)

   if (nargin < 2 || nargin > 3)
@@ -323,12 +298,31 @@

 endfunction

-%% test/octave.test/poly/residue-1.m
 %!test
 %! b = [1, 1, 1];
 %! a = [1, -5, 8, -4];
 %! [r, p, k, e] = residue (b, a);
-%! assert((abs (r - [-2; 7; 3]) < 1e-6
+%! assert ((abs (r - [-2; 7; 3]) < 1e-5
 %! && abs (p - [2; 2; 1]) < 1e-7
 %! && isempty (k)
 %! && e == [1; 2; 1]));
+%! k = [1 0];
+%! [b, a] = residue (r, p, k);
+%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
+%! && abs (a - [1, -5, 8, -4]) < 1e-12));
+
+%!test
+%! b = [1, 0, 1];
+%! a = [1, 0, 18, 0, 81];
+%! [r, p, k, e] = residue(b, a);
+%! r1 = [-5i; 12; +5i; 12]/54;
+%! r2 = conj(r1);
+%! p1 = [+3i; +3i; -3i; -3i];
+%! p2 = conj(p1);
+%! assert ((((abs (r - r1) < 1e-7) && (abs (p - p1) < 1e-7))
+%! ||       ((abs (r - r2) < 1e-7) && (abs (p - p2) < 1e-7)))
+%! && isempty (k)
+%! && (e == [1; 2; 1; 2]));
+%! [br, ar] = residue (r, p, k);
+%! assert ((abs (br - b) < 1e-12
+%! && abs (ar - a) < 1e-12));