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.
Copyright © 2006, Pascal Getreuer
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
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.
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.
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
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
Divide the exponential over to the right-hand side, and multiply by -(log b)/α to find
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  is the closed-form expression for iterated exponentiation:
Find z such that z z z ... = 2.
Solution: Using the formula,
The following command evaluates Wk(1) for k = −4,... 4:
w = lambertw((-4:4).',1)
These values are all solutions of w ew = 1. It is easy to verify numerically that they are solutions:
Problem 3 in the previous section mentioned the formula for evaluating iterated exponentiation:
For example, set z = 1.3, then its iterated exponentiation is approximately zlim = 1.4710.
zlim = -lambertw(-log(z))/log(z)
To verify this, computing z^z^z^...^z through 40 iterations shows that the iterated exponentiation does indeed converge to zlim.
zz(1) = z;
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
Largest error for b = -4: 2.51e-014
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 , 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.