bugGNU Octave - Bugs: bug #63937, mpoles uses relative, not...


bug #63937: mpoles uses relative, not absolute, tolerance for zero poles

Submitter:  Michael Asam <stromi>
Submitted:  Fri 17 Mar 2023 09:29:33 AM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  rik5
Originator Name:  Open/Closed:  * Closed
Release:  * 8.1.0 Operating System:  * Any
Fixed Release:  8.2.0 Planned Release:  8.2.0
* Mandatory Fields

Add a New Comment Rich Markup

Wed 22 Mar 2023 12:40:56 AM UTC, comment #5: 

I checked in the change on the stable branch of Octave (http://hg.savannah.gnu.org/hgweb/octave/rev/8eeb49ec16d8) which will become version 8.2 at the next bugfix release.  In the mean time, just replace mpoles.m in your local Octave installation with the file from this bug report.

Rik <rik5>
Group administrator
Tue 21 Mar 2023 03:45:44 PM UTC, comment #4: 

Hi Rik, your new mpoles() seems to work fine on my side. With your fix also residue() gives correct results in my case.

Many thanks!

Michael Asam <stromi>
Tue 21 Mar 2023 03:02:48 PM UTC, comment #3: 

I ended up re-writing most of the function.  Besides the change from relative to absolute tolerance for zero poles, I also added input validation (it was possible to enter a negative value for tolerance which would completely mess things up), BIST tests that verify compatibility with Matlab even for corner cases, BIST tests for input validation, re-wrote the documentation, changed variable names for clarity, and fixed incorrect indentation.  I've attached the new mpoles.m file for testing.

The only function that uses mpoles() in core is residue() and that function passes all of its BIST tests.  I don't know if any Octave packages use mpoles() and need to be tested.

(file #54518)

Rik <rik5>
Group administrator
Mon 20 Mar 2023 07:53:54 PM UTC, comment #2: 

This is actually two bugs in one.  The mpoles function needs to be fixed, but even with that in place residue also requires changes for the motivating example in this bug report.

I have modified mpoles.m, but while inspecting the code I found that it is very old (original check-in in 2007) and there are many improvements that should be made including using correct indentation, validating the inputs, having complete BIST tests, matching the variable names in the code to the documentation, and expanding the documentation to clearly state what the function does.  I've started on that.

Rik <rik5>
Group administrator
Fri 17 Mar 2023 03:41:51 PM UTC, comment #1: 

Your analysis looks accurate.

In your test, shouldn't the last input parameter to mpoles be 1 (true)?

If you have a patch for mpoles, I'm happy to review it and push it under your name.

Arun Giridhar <arungiridhar>
Group Member
Fri 17 Mar 2023 09:29:33 AM UTC, original submission:  

Here is an example where residue returns a wrong result if at least one pole is zero:

p = [0 -1e8 -9e6 -8e5 -8e3];
r = [10 -1 -2 -3 -4];
[b, a] = residue(r, p, []);

In this case I get the wrong result
b =
   4.0000e+00   5.5963e+08   7.0093e+15   7.2560e+21   5.7597e+25

The correct result for b would be:
b =
            0   1.2043e+08   3.0605e+15   4.3760e+21   5.7600e+25

The result for 'a' is okay.

After a closer look I found that the problem is within the function 'mpoles':
If a pole is nonzero, the function regards it unique, if the others are separated by a factor larger than 'tol' (same as Matlab does). 
But if a pole is zero, it regards it unique, if the other poles have a magnitude larger than 'tol' times the mean of all poles.

The latter is actually the root cause for the problem. Matlab in contrast treats 'tol' as an absolute tolerance for poles equal zero.

Here is an example for 'mpoles':

p = [0 -1e8 -9e6 -8e5 -8e3];
tol = 1e-3; % Hard coded in 'residue'
[mults, i] = mpoles(p, tol, 0)

The output is:
mults = 1  2  1  1  1
i     = 1  5  2  3  4

The output from Matlab is:
mults = 1  1  1  1  1
i     = 1  2  3  4  5

Can you change the behavior of 'mpoles', so that is treats 'tol' as an absolute tolerance for poles equal zero, like Matlab does?

Michael Asam <stromi>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #54518:  mpoles.m added by rik5 (5KiB - text/x-matlab)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by stromi (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.


    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-03-22 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
        Fixed ReleaseNone 8.2.0
        Planned ReleaseNone 8.2.0
    2023-03-21 rik5 Attached File- Added mpoles.m, #54518
        StatusIn Progress Patch Submitted
        Summaryresidue() returns incorrect result in case of pole at zero mpoles uses relative, not absolute, tolerance for zero poles
    2023-03-20 rik5 Summary`residue` returns incorrect result in case of pole at zero residue() returns incorrect result in case of pole at zero
    2023-03-20 rik5 StatusNone In Progress
        Assigned toNone rik5

    Back to the top

    Powered by Savane 3.12