[HP41] Lambert function RPN; question

[Albert Chan]

(12-27-2022 02:31 AM)Thomas Klemm Wrote:  It took me a while to figure out that I had used Halley's method with:

\(
\begin{align}
f(x) &= x e^x - a \\
f'(x) &= (1 + x) e^x \\
f''(x) &= (2 + x) e^x \\
\end{align}
\)

There was a problem with this setup. Derivatives may be huge, possibly bigger than f itself!
Had we use simple Newton's method, it may not converge at all.

>>> from mpmath import *
>>> a = 1e99
>>> x = ln(a) # W(a) guess ≈ 228
>>> def halley(x): y=exp(x); z=x*y; f=z-a; return f/(y+z - f/(y+z)*(y+z/2))
...
>>> for i in range(6): x -= halley(x); print x
...
225.982091142775
224.113871381155
222.808520559707
222.552199206554
222.550768955996
222.55076895575

Solve x*exp(x) = a, Halley's method, take 6 iterations.
Solve y*ln(y) = a, y=e^x, Newton's method only take 4.

>>> y = a # e^W(a) guess
>>> for i in range(4): y=(y+a)/(log(y)+1); print a/y
...
114.477962103205
222.273911298397
222.550768184143
222.55076895575

A better Halley setup is first temper down huge derivatives.
Bonus: Halley's correction term is simpler.

f = x + ln(x/a)
f' = 1 + 1/x
f'' = -1/x^2

>>> x = ln(a) # same guess
>>> def halley(x): f=x+ln(x/a); return f*x/(x+1 + f/(x+x+2))
...
>>> for i in range(2): x -= halley(x); print x
...
222.550764466154
222.55076895575