Numeric root-finding in Python

[Steven D'Aprano]

On Sun, 12 Feb 2012 12:18:15 -0800, Mark Dickinson wrote:

> On Feb 12, 6:41 am, Steven D'Aprano <steve
> +comp.lang.pyt... at pearwood.info> wrote:
> 
>>             err = -a/b  # Estimate of the error in the current w.
>>             if abs(err) <= 1e-16:
>>                 break
> 
> If the result you're expecting is around -1.005, this exit condition is
> rather optimistic:  the difference between the two Python floats either
> side of this value is already 2.22e-16, so you're asking for less than
> half a ulp of error!

I was gradually coming to the conclusion on my own that I was being 
overly optimistic with my error condition, although I couldn't put it 
into words *why*. Thanks for this Mark, this is exactly the sort of thing 
I need to learn -- as is obvious, I'm no expert on numeric programming.


> As to the rest;  your error estimate simply doesn't have enough
> precision.  The main problem is in the computation of a, where you're
> subtracting two almost identical values.  The absolute error incurred in
> computing w*exp(w) is of the same order of magnitude as the difference
> 'w*exp(w) - x' itself, so err has lost essentially all of its
> significant bits, and is at best only a crude indicator of the size of
> the error.  The solution would be to compute the quantities 'exp(w),
> w*exp(w), and w*exp(w) - x' all with extended precision.

Other than using Decimal, there's no way to do that in pure Python, is 
there? We have floats (double) and that's it.


> For the other
> quantities, there shouldn't be any major issues---after all, you only
> need a few significant bits of 'delta' for it to be useful, but with the
> subtraction that generates a, you don't even get those few significant
> bits.
> 
>> (The correct value for w is approximately -1.00572223991.)
> 
> Are you sure?  Wolfram Alpha gives me the following value for W(-1,
> -0.36787344117144249455719773322925902903079986572265625):
> 
> -1.005722239691522978...


I did say *approximately*. The figure I quote comes from my HP-48GX, and 
seems to be accurate to the precision offered by the HP.



-- 
Steven