The MATLAB function lambertw included in this package is a selfcontained Mfunction 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 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. Function Usagew = lambertw(z) computes the principal value of the Lambert W function, W_{0}. 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. BackgroundThe Lambert W function is defined as the function W(z) such that W(z) e^{W(z)} = z for all complex values z. As log z is the inverse of e^{z}, W(z) is the inverse of z e^{z}. 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 W_{k} for the kth branch. The principal branch, W_{0}(z), is realvalued for −1/e ≤ z. If −1/e ≤ z < 0, then the branch W_{−1}(z) is also realvalued. In the complex plane, a surface plot of W_{0}(z) is
Algebra with the Lambert W FunctionAlthough 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 buildingblock 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)e^{f(x)}. Problem 1. Solve y = (x − 1) e^{2x} for x. Solution: The righthand side is close to the necessary form f(x)e^{f(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) e^{2x}. This form is closer yet, the exponent is only missing the −2. Multiplying both sides by e^{−2} achieves the desired form, 2e^{−2}y = (2x − 2) e^{2x−2}. At this point, the Lambert W function can be applied to yield W(2e^{−2}y) = 2x − 2, which leads to the solution x = W(2e^{2}y)/2 + 1. Note that since W is multivalued, the solution is multivalued; there are multiple values of x satisfying the equation. Problem 2. Solve b^{x} = x^{α} for x. Solution: In this problem, there is initially no visible instance of the exponential function. Rewrite b^{x} as e^{x log b}, revealing
Divide the exponential over to the righthand side, and multiply by (log b)/α to find
The righthand side is now in the form f(x)e^{f(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 closedform expression for iterated exponentiation:
Find z such that z^{ z z ...} = 2. Solution: Using the formula,
DemosThe following command evaluates W_{k}(1) for k = −4,... 4: w = lambertw((4:4).',1) w = These values are all solutions of w e^{w} = 1. It is easy to verify numerically that they are solutions: w.*exp(w) ans = 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) zlim = To verify this, computing z^z^z^...^z through 40 iterations shows that the iterated exponentiation does indeed converge to zlim. zz(1) = z;
TestIdeally, 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.51e014 ImplementationThe Lambert W function is implemented numerically with approximations from series expansions followed by rootfinding. 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 rootfinder. As developed in [2], lambertw uses Halley's method, a fourthorder extension of Newton's rootfinding method. Convergence is very fast, usually requiring fewer than 5 iterations to reach machine accuracy. References

Research >