# best approximation to the LambertW(x) or exp(LambertW(x))

what is the best available approximation ( say up to 10 digits ) for LambertW(x) or exp(LambertW(x)) for x > 2000

-
Have you looked at functions.wolfram.com/ElementaryFunctions/ProductLog, for example? If so, can you be more explicit about what you want to know that isn't there? – Neil Strickland Mar 8 '11 at 14:28
Do you want an approximation in term of more elementary functions (then ten digits seems to be rather a lot) or just a way to calculate it numerically (then a root-finder seems to be the method of choice). – Fabian Mar 8 '11 at 16:34
Why don't you look at how GSL computes that function. – Zsbán Ambrus Mar 8 '11 at 16:37
– Shai Covo Mar 16 '11 at 14:44

The approximation described below is original, explicit (in some sense), and very accurate. It is closely related to this question (second paragraph).

So, you want to approximate the solution $W(x)$ of $x=W(x)e^{W(x)}$, for large $x$ (the order of $x$ does not play a very significant role in what follows). First, define $$\varphi (x,r) = 1 + \sum\limits_{k = 1}^{\left\lceil r \right\rceil } {\frac{{x^k [r - (k - 1)]^k }}{{k!}}} .$$ Now, consider the following series of approximations, where $r$ is assumed sufficiently large. The first one is $$\tilde W^1 (x,r) = \frac{1}{r}\ln \varphi (x,r).$$ Subsequent approximations are defined recursively by $$\tilde W^{n + 1} (x,r) = \frac{1}{r}\ln \bigg[\frac{{\tilde W^n (1 + \tilde W^n )}}{x}\varphi (x,r)\bigg].$$

Example. For $x=2000$, even $r$ as low as $80$ gives quite accurate results: $$\tilde W^5 (2000,80) \approx 5.83673149492073$$ and $$\tilde W^6 (2000,80) \approx 5.836731494908671,$$ while the exact solution (according to Wims Function Calculator) is $$W(2000) = 5.836731494908178747....$$ Based on many numerical results, this approximation seems quite interesting. Here is one further example. The Omega constant $\Omega$ is the value of $W(1)$: $$\Omega = W(1) \approx 0.5671432904097838729999686622.$$ With $r$ as low as $30$, we already get the following impressive approximations: $$\tilde W^1 (1,30) \approx 0.5710729200334063,$$ $$\tilde W^2 (1,30) \approx 0.5674569334624368,$$ $$\tilde W^3 (1,30) \approx 0.5671683899602143,$$ $$\tilde W^4 (1,30) \approx 0.5671452994467842,$$ $$\tilde W^5 (1,30) \approx 0.5671434512213455,$$ $$\tilde W^6 (1,30) \approx 0.5671433032818183,$$ $$\tilde W^7 (1,30) \approx 0.5671432914401158,$$ $$\tilde W^8 (1,30) \approx 0.567143290492256,$$ $$\tilde W^9 (1,30) \approx 0.5671432904163853,$$ $$\tilde W^{10} (1,30) \approx 0.5671432904103123,$$ $$\tilde W^{11} (1,30) \approx 0.5671432904098261,$$ $$\tilde W^{12} (1,30) \approx 0.5671432904097873.$$ So, $\tilde W^{12} (1,30) - W(1) \approx 3 \times 10^{-15}$. It is interesting to compare this sophisticated approximation with the standard one obtained from the converging sequence $\Omega_n \to \Omega$ defined by $\Omega_{n+1} = e^{-\Omega_n}$ (with initial value $\Omega_0$). For example, with $\Omega_0 = 0.5$, we only get $$\Omega_1 \approx 0.6065306597126334,$$ $$\Omega_2 \approx 0.545239211892605,$$ $$\Omega_3 \approx 0.5797030948780683,$$ $$\Omega_4 \approx 0.5600646279389019,$$ $$\Omega_5 \approx 0.5711721489772151,$$ $$\Omega_6 \approx 0.5648629469803235,$$ $$\Omega_7 \approx 0.5684380475700662,$$ $$\Omega_8 \approx 0.5664094527469208,$$ $$\Omega_9 \approx 0.5675596342622424,$$ $$\Omega_{10} \approx 0.5669072129354714,$$ $$\Omega_{11} \approx 0.5672771959707785,$$ $$\Omega_{12} \approx 0.5670673518537281.$$

-
Thaks a lot Covo. Please donot remove this answer. – Truth Seeker Mar 11 '11 at 13:29
– Shai Covo Mar 16 '11 at 14:43

Your question is not really clear about what you mean by 'approximation.' The Lambert W(x) function is implemented in various software packages, as ProductLog[x] in Mathematica, for example. And Mathematica can compute numerical values for specific x out to as many digits as you like:

N[ProductLog[2000],25]=5.836731494908178747954545

You may instead be asking about the asymptotic expansion of W(x) as x goes to $\infty$. Here the Mathematica input

FullSimplify[Normal[Series[ProductLog[x], {x, Infinity, 0}]], Assumptions -> x > E]

returns

$\log (x)-\log (\log (x))+\log (\log (x))/\log (x)+(\log (\log (x))-2) \log (\log (x))/(2 \log ^2(x))$.

So one may say that

$W(x) = \log (x)-\log (\log (x))+\log (\log (x))/\log (x)+O\left(1/\log(x)\right)$.

(If you're unfamiliar with asymptotic analysis and the Big Oh notation, you might start with http://en.wikipedia.org/wiki/Big_O_notation )

-
Mathematica is wrong here; the error is not O(1/x) but Theta((log log x/log x)^3). – Charles Mar 16 '11 at 14:38
@Charles: Thanks! I tried to edit to fix this. – Stopple Mar 23 '11 at 22:53

You want to approximate the solution $w:=W(x)$ of the equation $we^w =x$, respectively, the solution $u:=\exp(W(x))$ of of the equation $u\log(u)=x$, in dependence on the parameter $x\ge 2000$. In general, solutions of such equations are easily approximated by means of iterative methods. You may use the Newton method (check the linked article for the quadratic bounds on the approximation). Anyway, the Newton Method is very fast once you are conveniently close to the solution: I guess that the best thing here is to start with a more rough but more stable method, and pass to the Newton iteration as soon as the error became sufficiently small. Anyway, for these particular equations, of course, you may find ready the study on various approximations in all details: check e.g. the wiki article on the Lambert W function and the external links therein.

-

If you want a good approximation on a given interval, you can sometimes do significantly better than asymptotics using Chebyshev/Pade approximations or the Remez algorithm.

-
Hi Professor Israel. Thank you for this answer. You have been very helpful on the mapleprimes. I am so glad to see you here. – Truth Seeker Mar 25 '11 at 21:46

As Pietro says, $u=\exp(W(x))$ satisfies $u=x/\ln u$. This is a contraction mapping for large enough $u$, so just start with any old approximation, like $u=x/\ln x$, and do $u:=x/\ln u$ until it converges. For $x>1000$ it doesn't take more than about a dozen iterations to get 10 digits.

-

I don't have enough reputation to post comments. It appears that LambertW(x*(complex number)) where x=1,2,3,4,5,6,... can be calculated using tetration as in this Mathematica program:

Clear[nn, t, n, k, i];
nn = 75;
t[1, 1] = 1;
t[n_, k_] :=
t[n, k] =
If[n >= k,
Exp[-Sum[N[(1/ZetaZero[1])*t[n - i, k]], {i, 1, k - 1}]], (1/
ZetaZero[1])];
Table[t[nn, k]*(k - 1)*(1/ZetaZero[1]), {k, 1 + 1, 6 + 1}]


{0.00736674 - 0.0697969 I, 0.0235635 - 0.136089 I, 0.0466203 - 0.197008 I, 0.07436 - 0.251797 I, 0.104914 - 0.300521 I, 0.136894 - 0.343688 I}

Compared to Mathematicas ProductLog command:

Table[N[ProductLog[n*(1/ZetaZero[1])]], {n, 1, 6}]


{0.00736674 - 0.0697969 I, 0.0235635 - 0.136089 I, 0.0466203 - 0.197008 I, 0.07436 - 0.251797 I, 0.104914 - 0.300521 I, 0.136894 - 0.343688 I}

Link to question on mathematics stackexchange, LambertW(k)/k by tetration for natural numbers

I now noticed I do have enough reputation to post comments.

-