Re: Patch to residue.m

[Ben Abbott]

Unfortunately, I can't receive email from my usual server today :- 
( ... but did see your response on the octave web site.

So I thought I'd cut and paste a bit and respond.

> > Regarding, "toler" and "small" these two have two different  
> purposes.
> > "toler" is used to determine the pole multiplicity and "small" is  
> used
> > to round off to zero, pure real, or pure imaginary values. If  
> "toler"
> > is replace by "small" there will be no multiplicity of poles. So I
> > don't think that is practical.
>
> Well, that doesn't mean that small = toler = 1e-12.  Perhaps some  
> other value.  I'm just trying to justify in my mind the need for two  
> tolerance parameters vs. one.

I'd love to have a consistent and proper method to determine  
tolerances for determining multiplicity and zero/real/imag values. Do  
you have any idea how that might be done, or where to look for a  
mathematical method?

I spent a few hours googling yesterday in the hopes of finding a  
reference for who the error in the method used by roots.m might be  
estimated. Unfortunately, I did not find anything useful.

Today I decided to try some different examples to see how the errors  
vary with the order of the polynomial. I tried four cases.

(1) I assumed a varying order where all roots = pi.
(2) Same as (1), but when solving for the roots, I took the mean.
(3) I selected roots spiraling outward from the unit circle in the  
complex plane.
(4) I selected roots as p = 1:order

From these roots, I constructed representative polynomials, solved  
for the roots, and calculated the error in the solution.

I've attached the resulting plot.

While the plot doesn't offer much insight into how we might estimate  
the error of the roots, a few observations can be made.

(a) taking the mean to determine the pole value in the presence of  
multiplicity is a *good* idea!
(b) order and differences in values all effect the accuracy by which  
the roots are determined.

> > What we really need is numerically more accurate method for   
> determining
> > the roots. This entire discussion is the consequence of  inaccurate
> > roots. So if we are to address the ills of residue.m, a  better  
> roots.m
> > will be needed, in my opinion anyway.
>
> I think you are right.  The roots being off by o(1e-5) seems rather  
> poor.
> Averaging gets to o(1e-15) and machine precision is o(1e-16)?  One  
> would
> think maybe o(1e-10) or better is achievable.  Maybe you are trying  
> to improve
> on something that should be addressed elsewhere, and if that where  
> fixed your
> tolerance equations could be simplified.

Once an agreement is reached that averaging the pole values in the  
presence of multiplicity is proper, then there no longer remains is  
the errors introduced by roots.m ... well that any any adjustments we  
might want to apply in residue.m in consideration of those errors. By  
"adjustments", I refer to the conversion to zero/real/imag values  
based upon some tolerance.

There are more reliable methods to determine the roots of polynomials,  
but there are ***very*** slow :-(

If you're interested in looking over an example, Google "multroots".

Ben

roots_and_order.pdf
Description: Adobe PDF document