(28/48/50) Lambert W Function
|
03-27-2023, 08:57 PM
(This post was last modified: 03-27-2023 11:21 PM by Albert Chan.)
Post: #14
|
|||
|
|||
RE: (28/48/50) Lambert W Function
(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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)