(28/48/50) Lambert W Function

[Albert Chan]

(03-27-2023 04:45 PM)John Keith Wrote:  This line, h = x/(1+x) * s(x) returns zero or a very small number (usually 10^-12 or less)

Hi, John Keith

I assumed you meant correction around branch-point is small. That's normal.
My W guess formula is pretty accurate for a ≈ -1/e, both branches.

Solve s(x) = x + ln(x/a) = 0 (or its messy version) is probably a bad idea.
My previous post suggested it is better to solve for y = e^W(a), shown below.

lua> a = -0.3678 -- calculate W(a), 0 branch, by "hand"
lua> r, err = 0.36787944117144233, -1.2428753672788363E-17
lua> h = a + r + err
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.

lua> y = (y+a) / log1p((y-r-err)/r) -- Newton's step
lua> y, a/y, log(y)
0.37555110633147754    -0.9793607149578304      -0.9793607149578304

y's had doubled in accuracy, to full precision = e^W(a)

Quote:I have also seen several mentions of using Taylor series rather than iteration for values close to -1/e.

XCas> H := (1+x)*ln(1+x) - x
XCas> x0 := sqrt(2*H)               /* rough guess for above formula */
XCas> r0 := revert(series(x0,0,8))

\(\displaystyle x
+\frac{x^2}{6}
-\frac{x^3}{72}
+\frac{x^4}{270}
-\frac{23 x^5}{17280}
+\frac{19 x^6}{34020}
-\frac{11237 x^7}{43545600}
+\frac{13 x^8}{102060}\)

XCas> rs := series((r0/x)^2)      /* correction, inside square root */

\(\displaystyle
1+\frac{x}{3}
+\frac{x^3}{360}
-\frac{x^4}{810}
+\frac{23 x^5}{40320}
+ O\!\left(x^6\right)\)

This explained why my W guess formula is good (no x^2 terms, others tiny, with opposite sign)

Let's numerically check if series is good, with pade approximation matching to x^4
(replace 57 by 4032/71 ≈ 56.8, pade matched to x^5, but make no difference here)

lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3 + x^3/(360+160*x*(1-x/57)))
lua> y = r + r*x -- = e^W(a), if x is correct
lua> y, a/y, log(y)
0.3755511063314775      -0.9793607149578305      -0.9793607149578305

For comparison, my W code use simple guess (no bolded mess) + 1 Newton step.
Newton's step is very cheap; we had to convert x to W(a) guess anyway.

W(a) = a/y = a * (ln(y)+1)/(y+a), with y = r+r*x, we have:

Code:

x = a * log1p(x)/(r*x+h) -- Newton step

lua> W(-0.3678, nil, true) -- h=0 on the first try.
-0.9793607149578305
-0.9793607149578305