Numeric root-finding in Python
Steven D'Aprano
steve+comp.lang.python at pearwood.info
Sun Feb 12 18:05:47 EST 2012
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
More information about the Python-list
mailing list