Re: Patch to residue.m

[Ben Abbott]

On Monday, December 31, 2007, at 12:47PM, "Daniel J Sebald" <address@hidden> 
wrote:
>Doug Stewart wrote:
>
>>> 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.
>>>
>>> And considering the cast to zero after multiplicity is determined 
>>> seems worth investigating.
>>>
>>> Dan
>>>
>> For a quadratic there is an exact answer, and for cubics and for 
>> quartics:    but there is no exact method for 5th and higher polynomials.
>> Using residue one should recognize that it is a numerical approach to 
>> find approximate roots and residues. And if you recognize that they are 
>> approximate then the values used for tol and small seem appropriate.
>
>But why two?  And why of such drastic different values?

By two, you refer to "small" and to "toler" correct?

The reason there are two values, is because they have different purposes.

"small" is used to decide if a root's value could be set to zero/real/imag ... 
it is largely superficial. We could drop that part and the result would have 
little/no numerical consequence. However, it was present when I began working 
on this and others have expressed a desire that it remain ... so Iv'e left it 
alone.

"toler" is used to determine multiplicity, and needs to be relatively large due 
to the numerical instability/noise associated with poles of multiplicity. 

However, I'd be totally thrilled if someone could show me a improved 
mathematical method for root finding and/or determining multiplicity.

> 
>> Yes, if you are looking for roots in the 10e-12 range  then you you will 
>> have problems  - but this was true before the improvements.
>> 
>> If you want exact answers try Maxima, I did and found that it didn't 
>> even know how to do the exact 3 and 4th order answers..
>
>Not exact, but better, or at least consistent amongst methods.  It appears 
>that mpoles() has some code similar to the averaging technique.  Can that code 
>be reused somehow so that numerical effects are consistent?

I've looked, but have not come up with anything. For what its worth, I did 
compare the root finding in octave to matlab, and found them to be comparable. 
Actually, octave did a bit better; I assume because my build uses lapack/blas 
taylored for my architecture.

>
>Also, why is it important to have a tighter restriction on poles being zero?  
>If one simply uses the averaging technique and pole multiplicity, then there 
>might be, say, a multiplicity of 2 at 0.00000313 or something.

The tolerance for setting to zero/real/imag is not reliable. In *most* cases 
roots that should be pure real/imag or zero will be missed because the errors 
in their solution are much larger than "small".

>
>> If someone would write a better roots.m file we will all benefit. But  
>> it is not an easy task. ( I have been thinking about it since 
>> "Numerical  method" course in university in 1982) Using VPA (Variable 
>> Precision Arithmetic) yes you can get the roots to the nth decimal point 
>> but at a cost of cpu time-- so the conclusion that I have come to is:    
>> that this is a step in the right direction and is an improvement over 
>> what was here before.
>
>Sure, but so long as it is fresh in the mind, it'd be good to identify if and 
>where there is an issue if that means just leaving a comment in the code "root 
>accuracy needs to be improved".

Sure. Please ask to have it added to a list ... better yet you might try to 
take that on?

However, no need to add it to the code for residue. As residue isn't direclty 
impacted.

Ben