[Octave-bug-tracker] [bug #63937] `residue` returns incorrect result in case of pole at zero

[Michael Asam]

URL:
  <https://savannah.gnu.org/bugs/?63937>

                 Summary: `residue` returns incorrect result in case of pole
at zero
                   Group: GNU Octave
               Submitter: 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: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
                 Release: 8.1.0
         Discussion Lock: Any
        Operating System: Any
           Fixed Release: None
         Planned Release: None


    _______________________________________________________

Follow-up Comments:


-------------------------------------------------------
Date: Fri 17 Mar 2023 09:29:33 AM UTC By: Michael Asam <stromi>
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?







    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?63937>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/