Re: modified residues() for matlab compatibility

[John W. Eaton]

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