(28/48/50) Lambert W Function

[Albert Chan]

(03-27-2023 08:57 PM)Albert Chan Wrote:  lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3) -- better guess
lua> y = r + r*x -- = e^W(a), if x is correct
lua> y, a/y, log(y)
0.3755511062373703      -0.9793607152032429      -0.9793607152084147

W(a) 2 ways matched does not imply they are correct, only that we are close.
Interestingly, true W(a) may not be a "mean" of the two ways.

I was curious why the 2 W's do not bracket true W(a)
Doing a bit of error analysis, everything becomes clear.

Ignoring error of division, we have: relerr(x = a/y) ≈ - relerr(y)

ln(y*(1+ε)) ≈ ln(y) + ε
ln(y*(1+ε)) / ln(y) ≈ 1 + ε/ln(y)     relerr(ln(y)) ≈ relerr(y) / ln(y)

2 W's don't bracket true W(a) is because ln(y) ≈ x is negative.
2 W's have errors of the same sign.

relerr(x) + x*relerr(ln(y)) ≈ -relerr(y) + relerr(y) ≈ 0

lua> x, lny = -0.9793607152032429, -0.9793607152084147
lua> (x + x*lny) / (1+x)
-0.9793607149578298

Simply extrapolate for 0 error, we get good W(-0.3678) estimate.
Actually, extrapolate for 0 error is equivalent to a Newton's step!