pracma (version 1.9.9)

lambertWp: Lambert's W Function

Description

Principal real branch of the Lambert W function.

Usage

lambertWp(x)
lambertWn(x)

Arguments

x

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

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

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<="" code="">. 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.

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. http://www.apmaths.uwo.ca/~djeffrey/Offprints/W-adv-cm.pdf.

See Also

halley

Examples

##  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 * 1i