Lambert's W Function

Principal real branch of the Lambert W function.

lambertWp(x)
lambertWn(x)

Arguments

x

Numeric vector of real numbers >= -1/e.

Details

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.

Value

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).

References

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.

Note

See the examples how values for the second branch or the complex Lambert W function could be calculated by Newton's method.

See also

halley

Examples

##  Examples
lambertWp(0)          #=> 0
#> [1] 0
lambertWp(1)          #=> 0.5671432904097838...  Omega constant
#> [1] 0.5671433
lambertWp(exp(1))     #=> 1
#> [1] 1
lambertWp(-log(2)/2)  #=> -log(2)
#> [1] -0.6931472

# 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...
#> [1] 1.19183

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
# [1] -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05   NaN

if (FALSE) { # \dontrun{
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")
  } # }

## 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
# [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
#> $zero
#> [1] -3.577152e+00  9.544935e-18
#> 
#> $fnorm
#> [1] 1.389481e-17
#> 
#> $niter
#> [1] 6
#> 
newtonsys(F, c(-1, -1), z0 = -pi/2)  #=> -1.5707963267949i = -pi/2 * 1i
#> $zero
#> [1]  4.436124e-18 -1.570796e+00
#> 
#> $fnorm
#> [1] 1.006197e-16
#> 
#> $niter
#> [1] 8
#>