Research‎ > ‎

Lambert W Function Implementation


Download lambertw

The MATLAB function lambertw included in this package is a self-contained M-function implementing the Lambert W function, also known as the Ω function or "product log." It is equivalent to the numerical mode of the Symbolic Math Toolbox's lambertw, however, no toolboxes are required to use this function.


License (BSD)

Copyright © 2006, Pascal Getreuer
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  • Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
  • Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Function Usage

w = lambertw(z) computes the principal value of the Lambert W function, W0. The input z may be a complex scalar or array. For real z, the result is real on the principal branch for z ≥ −1/e.

w = lambertw(b,z) specifies which branch of the Lambert W function to compute. If z is an array, b may be either an integer array of the same size as z or an integer scalar. If z is a scalar, b may be an array of any size.

Background

The Lambert W function is defined as the function W(z) such that

W(z) eW(z) = z

for all complex values z. As log z is the inverse of ez, W(z) is the inverse of z ez. Like the complex logarithm, the Lambert W function is multivalued with a countably infinite number of branches. The branches are enumerated by the integers and are conventionally denoted by Wk for the kth branch.

The principal branch, W0(z), is real-valued for −1/e ≤ z. If −1/e ≤ z < 0, then the branch W−1(z) is also real-valued. In the complex plane, a surface plot of |W0(z)| is

x = linspace(-6,6,51);
y = linspace(-6,6,51);
[x,y] = meshgrid(x,y);
w = lambertw(x + i*y);

surf(x,y,abs(w));
axis([-6,6,-6,6,0,2.5]);
view(40,32);
xlabel('Re z');
ylabel('Im z');
title('|W_0(z)|');

Algebra with the Lambert W Function

Although the Lambert W function may not be so widely known as the inverse trigonometry functions, it has essentially the same purpose as acos, asin, etc. as a building-block tool for solving equations. The general approach is to manipulate all occurrences of the unknown variable x into an expression of the form f(x)ef(x).

Problem 1. Solve y = (x − 1) e2x for x.

Solution: The right-hand side is close to the necessary form f(x)ef(x), but some manipulation is necessary to change the factor (x − 1) and the exponent 2x into the same expression. Multiplying both sides by 2,

2y = (2x − 2) e2x.

This form is closer yet, the exponent is only missing the −2. Multiplying both sides by e−2 achieves the desired form,

2e−2y = (2x − 2) e2x−2.

At this point, the Lambert W function can be applied to yield

W(2e−2y) = 2x − 2,

which leads to the solution x = W(2e-2y)/2 + 1. Note that since W is multivalued, the solution is multivalued; there are multiple values of x satisfying the equation.

Problem 2. Solve bx = xα for x.

Solution: In this problem, there is initially no visible instance of the exponential function. Rewrite bx as ex log b, revealing

ex log b = xα
ex (log b)/α = x.

Divide the exponential over to the right-hand side, and multiply by -(log b)/α to find

1 = x e-x (log b)/α
-(log b)/α = (-x (log b)/α) e-x (log b)/α.

The right-hand side is now in the form f(x)ef(x) where the Lambert W function can be applied:

W(-(log b)/α) = -x (log b)/α.

The solution is x = −(α / log b) W(−(log b)/α).

Problem 3. An old result [1] is the closed-form expression for iterated exponentiation:

z z z ... =   W(−log z) .
−log z

Find z such that z z z ... = 2.

Solution: Using the formula,

W(−log z) / log z = 2
W(−log z) = −2 log z
−log z = −2 log z e−2 log z
e2 log z = 2
z = 21/2.

Demos

The following command evaluates Wk(1) for k = −4,... 4:

w = lambertw((-4:4).',1)
w =

-3.1630 -23.4277i
-2.8536 -17.1135i
-2.4016 -10.7763i
-1.5339 - 4.3752i
0.5671
-1.5339 + 4.3752i
-2.4016 +10.7763i
-2.8536 +17.1135i
-3.1630 +23.4277i

These values are all solutions of w ew = 1. It is easy to verify numerically that they are solutions:

w.*exp(w)
ans =

1.0000 - 0.0000i
1.0000 - 0.0000i
1.0000 + 0.0000i
1.0000
1.0000
1.0000
1.0000 - 0.0000i
1.0000 + 0.0000i
1.0000 + 0.0000i

Problem 3 in the previous section mentioned the formula for evaluating iterated exponentiation:

z z z ... =   W(−log z) .
−log z

For example, set z = 1.3, then its iterated exponentiation is approximately zlim = 1.4710.

zlim = -lambertw(-log(z))/log(z)
zlim =

1.4710

To verify this, computing z^z^z^...^z through 40 iterations shows that the iterated exponentiation does indeed converge to zlim.

zz(1) = z;
zlim = -lambertw(-log(z))/log(z);

for k = 1:40
zz(k+1) = z^zz(k);
end

plot(zz,'.-');
  k    zz     zlim-zz
1 1.4065 0.065
2 1.4463 0.025
3 1.4615 0.0095
4 1.4673 0.0037
5 1.4696 0.0014
6 1.4704 0.00055
7 1.4708 0.00021
8 1.4709 8.1e-005
9 1.4710 3.1e-005
10 1.4710 1.2e-005
20 1.4710 8.9e-010
30 1.4710 6.5e-014
40 1.4710 2.2e-016
 

Test

Ideally, lambertw(b,z)*exp(lambertw(b,z)) = z for any complex z and any integer branch index b, but this is limited by machine precision. The inversion error |lambertw(b,z)*exp(lambertw(b,z)) - z| is small but worth minding.

Experimentation finds that the error is usually on the order of |z|×10−16 on the principal branch. This test computes the inversion error over the square [−10,10]×[−10,10] in the complex plane, large enough to characterize the error away from the branch points at z = 0 and −1/e.

N = 81;    % Use NxN points to sample the complex plane

R = 10; % Sample in the square [-R,R]x[-R,R]
x = linspace(-R,R,N);
y = linspace(-R,R,N);
[xx,yy] = meshgrid(x,y);
z = xx + 1i*yy;

for b = -4:4
w = lambertw(b,z);
InvError = abs(w.*exp(w) - z);
fprintf('Largest error for b = %2d: %.2e\n',b,max(InvError(:)));
end
Largest error for b = -4:  2.51e-014
Largest error for b = -3: 2.39e-014
Largest error for b = -2: 1.39e-014
Largest error for b = -1: 7.94e-015
Largest error for b = 0: 5.40e-015
Largest error for b = 1: 7.94e-015
Largest error for b = 2: 1.39e-014
Largest error for b = 3: 2.39e-014
Largest error for b = 4: 2.51e-014

Implementation

The Lambert W function is implemented numerically with approximations from series expansions followed by root-finding. Depending on the desired branch and the proximity to the branch points at z = 0 and −1/e, different series expansions are used as initializations to the root-finder.

As developed in [2], lambertw uses Halley's method, a fourth-order extension of Newton's root-finding method. Convergence is very fast, usually requiring fewer than 5 iterations to reach machine accuracy.

References

[1] G. Eisenstein, "Entwicklung von αα...." J. reine angewandte Math., vol. 28, 1844.
[2] R.M. Corless, G.H. Gonnet, D.E.G. Hare, G.J. Jeffery, and D.E. Knuth. "On the Lambert W Function." Advances in Computational Mathematics, vol. 5, 1996.