DEFINITION OF THE LAMBERT W FUNCTION
The Lambert W function is a rather unique and useful function that satisfies a special property.
For every non-negative real number $z$, there exists a unique non-negative real number $x$ such that $$x\exp(x) = z.$$
If this condition is satisfied, we say that $x$ is the Lambert W of $z$, and we denote it by $$W_0(z) = x.$$
This function, $W_0$, is the “principal branch” of the Lambert W function. It accepts non-negative real numbers as inputs and provides non-negative real numbers as outputs. We express this in mathematical notation as $W_0:[0,\infty) \rightarrow [0, \infty)$. Alternatively, $W_0$ can be defined as the inverse of the function $$g(x)=x\exp(x).$$
LAMBERT GROWTH RATE
We can use the Lambert W function to solve the delayed growth differential equation $$y'(t) = \beta y(t-T)$$ where $\beta$ and $T$ are positive real numbers.
It is not difficult to show that $$y(t) = \exp(\alpha T)$$ solves the delayed growth differential equation if and only if $$\alpha = W_0(\beta T)/T.$$
We define the Lambert Growth Rate for $\beta$ and $T$ to be $$\mathrm{Lambda\ Growth\ Rate}=\alpha = W_0(\beta T)/T.$$ (Equivalently, $\alpha$ is the unique real number such that the derivative of $y(t) =\exp(\alpha t)$ is $y(t-T)$.)
CALCULATION OR ESTIMATION
Often, we don’t have the Lambert W function readily available on our calculators or programming environments. So, how do we calculate or approximate it? Let’s look at some methods to compute or estimate this function. Let’s assume $\beta = 0.2$ and $T=30$ for the methods below.
The various methods for calculating or approximating the Lambert Growth Rate are as follows:
- Wolfram Alpha: This is a computational knowledge engine that can be used to solve the equation above directly. You could type
"solve a*exp(a*30) = 0.2"into the input box of Wolfram Alpha. The answer given by it would be approximately 0.0477468. Alternatively, to calculate $W_0( \beta T )/T$, you could type"W_0(0.2*30)/30". - Python and Scipy: Python is a popular programming language and scipy is a library for scientific computation in Python. You can use scipy’s implementation
from scipy.special import lambertw. - Excel: You can access the Lambert W function with the MoreFunc add-in for Excel.
- JavaScript: In Javascript, you can use the math library’s implementation of the Lambert W function. The code would be
math.lambertW(x). - Approximations: If you can accept an error of 1 or 2 percent, there are a few approximations available. For instance, if $1.6 \leq x \leq 22$, the Lambert W function can be approximated by the function $$w_1(x) = (.5391 x – .4479)^{(1/2.9)}.$$ We can use this approximation to estimate the growth rate
$$W_0( 30*0.2)/30= W_0(6)/30\approx (.5391\times 6 – .4479)^{(1/2.9)}/30\approx 0.047463321.$$ Also, if $0\leq x \leq 2$, another approximation function is $$w_2(x) = \frac{x}{(1+x)(1-0.109 x)}.$$ - Bisection Method: We want to find the value of $\alpha$ where $$\alpha \exp(\alpha T) = \beta.$$ That value of $\alpha$ is the “root” of the function $$f(\alpha)=\alpha \exp(\alpha T) – \beta.$$ A root of a function $f(\alpha)$ is any $\alpha$ where $f(\alpha)=0$. The bisection method continuously splits an interval in half to narrow down the root of a function. It assumes the function changes sign around the root, and iteratively refines the search interval. In our case, $f(0)= -\beta<0$ and $$f(\beta) = \beta\exp(T\beta) -\beta>0,$$ so we know that $\alpha=0$ is too low and $\alpha=\beta$ is too high. The bisection method would now test the midpoint of the interval $[0,\beta]$ which is $\beta/2.$ If $f(\beta/2)>0$ the bisection method begins again on the interval $[0,\beta/2]$ and if $f(\beta/2)<0$ the bisection method would continue with the interval $[\beta/2, \beta]$. Every iteration cuts the size of the interval in half.
- Newton-Raphson Method: The Newton-Raphson method is a popular root-finding algorithm that uses tangents to the function to find the roots. It can be used to improve the accuracy of the approximations above. The formula for estimating the solution of $f(x)=0$ is $$n(x) = x – f(x)/f'(x).$$For our problem $$n(x) = x-\frac{x e^{T x}-\beta}{(T x +1)e^{T x}}.$$
Simulation
In the game Master of Orion, the economic growth rate in the early part of the game can be estimated by $\alpha=W_0(T\frac{I}{C})$ where $I$ represents the income of a “mature” planet, $C$ is the cost of constructing a colony ship, $T$ is the number of years between the construction year of a colony ship and the year that the colonized planet is “mature”. I ran a simulation using typical but random values where $10\leq I\leq200$, $200\leq C\leq575$, and $10\leq T\leq90$. Every time I generated random $I$, $C$, and $T$ values, I would use either $w_1$ or $w_2$ to estimate the Lambert growth rate $W_0(TI/C)/T$. Let $x= TI/C$. If $x<1.8$, I used $w_2$ and I used $w_1$ otherwise. The average error using $w_1$ and $w_2$ was 0.0016. If I then applied $n$ to the approximation, then the average error was 0.00017, about 10 times more accurate. When I applied $n$ twice the average error dropped to 0.000005 and the worst case error was 0.0001. If you want 16 digits of accuracy, you need to apply the function $n$ five times.
Applying the Newton Raphson to estimate the Lambert growth rate for $T=30$ and $I/C = \beta = 0.2$ and, gives
$$
\begin{aligned}
W_0(I T/C)/C = W_0(6)/30 \approx w_1(6)/30 &\approx 0.047463321\\
W_0(6)/30 \approx n(w_1(6)/30) &\approx 0.047748535122\\
W_0(6)/30 \approx n(n(w_1(6)/30)) &\approx 0.047746825925 \\
W_0(6)/30&\approx 0.047746825863
\end{aligned}
$$
Summary
In this post we presented several methods for approximating the Lambert Growth Rate $$W_0(\beta T)/T.$$ If you don’t have Wolfram Alpha or access to a programming environment that includes the Lambert W function, then one of the best methods for finding the solution is the Bisection Method. If $x=\beta T<100$, then using $w_1$ or $w_2$ approximates the solution to within 2%. One iteration of Newton-Raphson typically reduces the error by a factor of 10. More iterations of Newton-Raphson significantly improve the approximation.