在 2012年2月12日星期日UTC+8下午2时41分20秒，Steven D'Aprano写道：
> This is only peripherally a Python problem, but in case anyone has any
> good ideas I'm going to ask it.
>
> I have a routine to calculate an approximation of Lambert's W function,
> and then apply a root-finding technique to improve the approximation.
> This mostly works well, but sometimes the root-finder gets stuck in a
> cycle.
>
> Here's my function:
>
> import math
> def improve(x, w, exp=math.exp):
> """Use Halley's method to improve an estimate of W(x) given
> an initial estimate w.
> """
> try:
> for i in range(36): # Max number of iterations.
> ew = exp(w)
> a = w*ew - x
W*EXP(W) can converge for negative values of W
> b = ew*(w + 1)
b=exp(W)*W+W
> err = -a/b # Estimate of the error in the current w.
What's X not expalained?
> if abs(err) <= 1e-16:
> break
> print '%d: w= %r err= %r' % (i, w, err)
> # Make a better estimate.
> c = (w + 2)*a/(2*w + 2)
> delta = a/(b - c)
> w -= delta
> else:
> raise RuntimeError('calculation failed to converge', err)
> except ZeroDivisionError:
> assert w == -1
> return w
>
>
> Here's an example where improve() converges very quickly:
>
> py> improve(-0.36, -1.222769842388856)
> 0: w= -1.222769842388856 err= -2.9158979924038895e-07
> 1: w= -1.2227701339785069 err= 8.4638038491998997e-16
> -1.222770133978506
>
> That's what I expect: convergence in only a few iterations.
>
> Here's an example where it gets stuck in a cycle, bouncing back and forth
> between two values:
>
> py> improve(-0.36787344117144249, -1.0057222396915309)
> 0: w= -1.0057222396915309 err= 2.6521238905750239e-14
> 1: w= -1.0057222396915044 err= -2.6521238905872001e-14
> 2: w= -1.0057222396915309 err= 2.6521238905750239e-14
> 3: w= -1.0057222396915044 err= -2.6521238905872001e-14
> 4: w= -1.0057222396915309 err= 2.6521238905750239e-14
> 5: w= -1.0057222396915044 err= -2.6521238905872001e-14
> [...]
> 32: w= -1.0057222396915309 err= 2.6521238905750239e-14
> 33: w= -1.0057222396915044 err= -2.6521238905872001e-14
> 34: w= -1.0057222396915309 err= 2.6521238905750239e-14
> 35: w= -1.0057222396915044 err= -2.6521238905872001e-14
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> File "<stdin>", line 19, in improve
> RuntimeError: ('calculation failed to converge', -2.6521238905872001e-14)
>
> (The correct value for w is approximately -1.00572223991.)
>
> I know that Newton's method is subject to cycles, but I haven't found any
> discussion about Halley's method and cycles, nor do I know what the best
> approach for breaking them would be. None of the papers on calculating
> the Lambert W function that I have found mentions this.
>
> Does anyone have any advice for solving this?
>
>
>
> --
> Steven
I sugest you can use Taylor's series expansion to speed up w*exp(w)
for negative values of w.
--
http://mail.python.org/mailman/listinfo/python-list
