The Lambert W function

YouTube’s algorithm kept presenting me with videos of equations whose solutions could be expressed in terms of the Lambet W function. The function, named after Johann Heinrich Lambert (1728 to 1777), is specified by:

W(z)e^{W(z)}=z, \quad z \in \mathbb{C}

W is multivalued, even when z and W are real. W(x) has a single value when x \ge 0, x \in \mathbb{R}. That is referred to as the principal branch and denoted W_0.

As:

0 \times e^0 = 0

W_0(0) = 0.

There are two values of W when -1/e \lt x \lt 0. One value is on the principal branch and W_0(x) \gt -1. The other value is on a branch denoted W_{-1} and W_{-1}(x) \lt -1.

As:

(-1)e^{-1}=-1/e

W(-1/e)=-1. That point is considered to be on both the principal branch and W_{-1}.

You can see that there are two values by identifying that there is a single maximum or minimum of x=We^W:

\frac {d}{dW} (We^W) = e^W+We^W=e^W(1+W)=0 \Rightarrow W=-1

and that:

\lim_{W \to -\infty} x = 0^-, \quad \lim_{W \to \infty} x = \infty

Haskell

A Hoogle search suggests that the hmatrix-special package provides lambert_W0 :: Double -> Double. That is a wrapper around a C function provided by the GNU Scientific Library (GSL).

The Flint2 package provides acb_lambertw :: Ptr CAcb -> Ptr CAcb -> Ptr CFmpz -> CInt -> CLong -> IO (). That is a wrapper around a C function provided by the FLINT: Fast Library for Number Theory.

Iacono, R., Boyd, J.P. New approximations to the principal real-valued branch of the Lambert W-function. Adv Comput Math 43, 1403–1436 (2017) provides an iteration to approximate W_0(x):

w_{n+1}(x) = \frac {w_n(x)} {1+w_n(x)} \left( 1+\ln \left( \frac {x} {w_n(x)}\right) \right), \quad xw_n(x) \gt 0, w_n(x) \neq -1

If w_n(x) = W_0(x), then:

w_{n+1}(x) = \frac {W_0(x)} {1+W_0(x)} \left( 1+\ln \left( \frac {x} {W_0(x)}\right) \right) = \frac {W_0(x)} {1+W_0(x)} \left( 1 + W_0(x) \right) = W_0(x)

Following Lóczi, Lajos. (2022). Guaranteed- and high-precision evaluation of the Lambert W function. Applied Mathematics and Computation. 433. 10.1016/j.amc.2022.127406, recommended inital values are:

w_0(x) = \begin{cases} \dfrac {ex \ln (1 + \sqrt{1+ex})} {ex + 1 + \sqrt{1+ex}} & -1/e < x < 0,\\ x/e & 0 < x < e,\\ \ln x - \ln {\ln x} & e < x \end{cases}

With those initial values:

w_{n+1}(x) \gt w_n(x), \quad x \gt 0 w_{n+1}(x) \lt w_n(x), \quad x \lt 0

For 0 < x < e, w_1(x) is:

\frac {x/e} {1 + x/e} \left( 1 + \ln \frac {x} {x/e} \right) = \frac {2x} {x + e}

As that is easy to calculate, it seems a better initial value.

This can be coded as a total function as follows:

Example

An example equation to solve from a YouTube video:

2^x + x = 5, \quad x \in \mathbb{R}

If 5 - x = y, then:

2^{5-y} = 2^5 2^{-y} = y

2^5 = 32 = y2^y = ye^{y \ln 2}

32 \ln 2 = y \ln 2 e^{y \ln 2}

So:

W_0(32 \ln 2) = y \ln 2

y = \frac {W_0(32 \ln 2)} {\ln 2}

and:

x = 5 - \frac {W_0(32 \ln 2)} {\ln 2}

An empirical test using stack ghci: