Numeric root-finding in Python

[Mark Dickinson]

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!

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.  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...

so it looks as though the values you're getting are at least
alternating around the exact value.


--
Mark