Numeric root-finding in Python
mdickinson at enthought.com
Sun Feb 12 21:18:15 CET 2012
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:
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,
so it looks as though the values you're getting are at least
alternating around the exact value.
More information about the Python-list