Lambert's W Function
Principal real branch of the Lambert W function.
lambertWp(x) lambertWn(x)
x |
Numeric vector of real numbers |
The Lambert W function is the inverse of x --> x e^x, with two
real branches, W0 for x >= -1/e and W-1 for -1/e <= x < 0.
Here the principal branch is called lambertWp, tho other one
lambertWp, computed for real x.
The value is calculated using an iteration that stems from applying Halley's method. This iteration is quite fast and accurate.
The functions is not really vectorized, but at least returns a vector of
values when presented with a numeric vector of length >= 2.
Returns the solution w of w*exp(w) = x for real x
with NaN if x < 1/exp(1) (resp. x >= 0 for the
second branch).
See the examples how values for the second branch or the complex Lambert W function could be calculated by Newton's method.
Corless, R. M., G. H.Gonnet, D. E. G Hare, D. J. Jeffrey, and D. E. Knuth (1996). On the Lambert W Function. Advances in Computational Mathematics, Vol. 5, pp. 329-359.
## Examples
lambertWp(0) #=> 0
lambertWp(1) #=> 0.5671432904097838... Omega constant
lambertWp(exp(1)) #=> 1
lambertWp(-log(2)/2) #=> -log(2)
# The solution of x * a^x = z is W(log(a)*z)/log(a)
# x * 123^(x-1) = 3
lambertWp(3*123*log(123))/log(123) #=> 1.19183018...
x <- seq(-0.35, 0.0, by=0.05)
w <- lambertWn(x)
w * exp(w) # max. error < 3e-16
# [1] -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 NaN
## Not run:
xs <- c(-1/exp(1), seq(-0.35, 6, by=0.05))
ys <- lambertWp(xs)
plot(xs, ys, type="l", col="darkred", lwd=2, ylim=c(-2,2),
main="Lambert W0 Function", xlab="", ylab="")
grid()
points(c(-1/exp(1), 0, 1, exp(1)), c(-1, 0, lambertWp(1), 1))
text(1.8, 0.5, "Omega constant")
## End(Not run)
## Analytic derivative of lambertWp (similar for lambertWn)
D_lambertWp <- function(x) {
xw <- lambertWp(x)
1 / (1+xw) / exp(xw)
}
D_lambertWp(c(-1/exp(1), 0, 1, exp(1)))
# [1] Inf 1.0000000 0.3618963 0.1839397
## Second branch resp. the complex function lambertWm()
F <- function(xy, z0) {
z <- xy[1] + xy[2]*1i
fz <- z * exp(z) - z0
return(c(Re(fz), Im(fz)))
}
newtonsys(F, c(-1, -1), z0 = -0.1) #=> -3.5771520639573
newtonsys(F, c(-1, -1), z0 = -pi/2) #=> -1.5707963267949i = -pi/2 * 1iPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.