A recurrence relation which approximates the Lambert W Function

Question

I have the following recurrence relation, which seems to approximate the Lambert W Function pretty well:

$W(x)\approx\log(f_n(x))$ for a sufficiently large value of $n$, where:

• $f_0(x) = x$
• $f_{k+1}(x) = x/\log(f_k)$

I have verified it using Python:

from mpmath import mp
from mpmath import lambertw
from decimal import Decimal
from decimal import getcontext

getcontext().prec = mp.dps = 80

def W(x):
    f = Decimal(x)
    for n in range(100):
        f = x / f.ln()
    return f.ln()

for x in range(3,10):
    print(lambertw(x))
    print(W(x))

The problem is that it converges to a real value only for $x\geq e$. And since the Lambert W Function has a known asymptotic expansion which converges for $-1/e\leq x\leq1/e$, I am essentially interested only in the range $1/e<x<e$.

Is there a way for me to tweak the above in order to support this range (for example, by choosing a different initial value)?

---

Though not necessarily required in this context, here is some further analysis of this recurrence relation:

• For $x\geq e$ is seems to converge well
• For $1<x<e$, it reaches a negative value at some point (when choosing a large enough value of $n$), and then turns complex (due to applying $\log$ on a negative value)
• For $x=1$, it yields division by $0$ of course
• For $0<x<1$, it reaches a negative value immediately, and then turns complex (due to applying $\log$ on a negative value)
• For $x=0$, it yields division by $-\infty$ of course
• For $x<0$, it turns complex immediately (due to applying $\log$ on a negative value)
• For $x<-1/e$, it is no longer relevant since $W(x)$ has no real solutions

Thank you all!!!

Answer 1

Your algorithm is basically based on the continued fraction expansion found in Wiki for the Lambert W function's principal branch. For this, consider the equivalent iteration (1): $g_{k+1}(x)=\log(x/g_k(x))$. Assuming it converges to a fixed point $x_0$, the iteration will satisfy:

$$x_0=\log\left(\frac{x}{x_0}\right)\Leftrightarrow$$ $$\exp(x_0)\cdot x_0=x\Leftrightarrow$$ $$x_0=W(x)$$

In order for the iterates of the map: $g(x)=\log\left(\frac{x}{g_0(x)}\right)$ to converge to the fixed point $x_0$, the multipler of the map $g$ has to be:

$$|g'(x_0)|=\left |\frac{1}{(x/g_0(x))}\cdot\frac{1}{g_0(x)}\right|_{x=x_0}=\frac{1}{|x_0|}=\frac{1}{|W(x)|}\lt 1$$

In other words you need: $|W(x)|\gt 1$ or equivalently $|x|\gt e$. Otherwise the fixed point $x_0$ might be repulsive. Combined with the domain of the principal branch of $W$ which is $[-1/e,\infty)$, you get $x\gt e$ (It may be that Wiki's continued fraction expansion has a typo ([2]): It says it's valid for $|W(x)|>e$, which gives $x>e^{e+1}$ which is a subset of where this iteration converges.

Concluding, if you want a similar algorithm for $x\in(-1/e,e)$, use the continued fraction expansion above the one quoted in Wiki, which is equivalent to the iteration:

$$f_{k+1}(x)=\frac{x}{\exp(f_k(x))}, k\in\mathbb{N}$$

which is valid for $|W(x)|<1\Leftrightarrow x\in(-1/e,e)$.

1: Unless I have a typo somewhere, note that the iteration $g_k$ is equivalent to your iteration $f_k$ and they are related as: $g_n(x)=\log(f_n)$, with $g_0(x)=\log(f_0(x))=\log(x)$.

[2]: It looks indeed as if it is a typo: Check the recursive expressions (73) on p.203 in the quoted source:

$\begin{align} v_{n+1}&=v_n+p_n\\ p_{n+1}&=-\log(1+p_n/v_n) \end{align}$

If you set: $v_0=\log(z)$ and $p_0=-\log(\log(z))$, then unwinding the recursion for $v_n$ gives the continued fraction expansion in Wiki. But in the source (right top of p. 203) it is mentioned that $v_n\to W(z)$ if $|W(z)|>1\Leftrightarrow |x|>e$.

Answer 2

You're essentially using the iteration

$$g_n\exp(g_{n+1})=x\implies g_{n+1}=\ln\left(\frac x{g_n}\right)$$

which converges for large $x$. Inverting this so that we have

$$g_{n+1}\exp(g_n)=x\implies g_n=x\exp(-g_n)$$

will converge for $-1/e\le x\le e$.

A good initial approximation is $g_0=\ln(x+1)$, which works well for for both small and large values.

Faster methods exist, such as Newton's method, which leads to

$$g_{n+1}=g_n-\frac{g_n\exp(g_n)-x}{(g_n+1)\exp(g_n)}=\frac{g_n^2+x\exp(-g_n)}{g_n+1}$$

or Aitken's acceleration method:

$$g_{n+1}=g_n-\frac{(g_n'-g_n)^2}{g_n''-2g_n'+g_n}\\g_n'=x\exp(-g_n)\\g_n''=x\exp(-g_n')$$

which both converge for all $x\ge-1/e$.