# [R] Lambert's W function

Ben Bolker bolker at zoo.ufl.edu
Wed Nov 26 15:52:44 CET 2003

  This implementation, originally written by Nici Schraudolph, allows you
to choose which branch you want.  I've checked the answers for complex
arguments, non-systematically, against Mathematica's ProductLog function.

Ben Bolker

lambertW = function(z,b=0,maxiter=10,eps=.Machine$double.eps, min.imag=1e-9) { if (any(round(Re(b)) != b)) stop("branch number for W must be an integer") if (!is.complex(z) && any(z<0)) z=as.complex(z) ## series expansion about -1/e ## ## p = (1 - 2*abs(b)).*sqrt(2*e*z + 2); ## w = (11/72)*p; ## w = (w - 1/3).*p; ## w = (w + 1).*p - 1 ## ## first-order version suffices: ## w = (1 - 2*abs(b))*sqrt(2*exp(1)*z + 2) - 1 ## asymptotic expansion at 0 and Inf ## v = log(z + as.numeric(z==0 & b==0)) + 2*pi*b*1i; v = v - log(v + as.numeric(v==0)) ## choose strategy for initial guess ## c = abs(z + exp(-1)); c = (c > 1.45 - 1.1*abs(b)); c = c | (b*Im(z) > 0) | (!Im(z) & (b == 1)) w = (1 - c)*w + c*v ## Halley iteration ## for (n in 1:maxiter) { p = exp(w) t = w*p - z f = (w != -1) t = f*t/(p*(w + f) - 0.5*(w + 2.0)*t/(w + f)) w = w - t if (abs(Re(t)) < (2.48*eps)*(1.0 + abs(Re(w))) && abs(Im(t)) < (2.48*eps)*(1.0 + abs(Im(w)))) break } if (n==maxiter) warning(paste("iteration limit (",maxiter, ") reached, result of W may be inaccurate",sep="")) if (all(Im(w)<min.imag)) w = as.numeric(w) return(w) } --------- \name{lambertW} \alias{lambertW} \title{Lambert W function} \description{ Computes the Lambert W function, giving efficient solutions to the equation x*exp(x)==x } \usage{ lambertW(z, b = 0, maxiter = 10, eps = .Machine$double.eps, min.imag = 1e-09)
}
\arguments{
\item{z}{(complex) vector of values for which to compute the function}
\item{b}{(integer) vector of branches: b=0 specifies the principal
branch, 0 and -1 are the ones that can take non-complex values}
\item{maxiter}{maximum numbers of iterations for convergence}
\item{eps}{convergence tolerance}
\item{min.imag}{maximum magnitude of imaginary part to chop when
returning solutions}
}
\details{
Compute the Lambert W function of z.  This function satisfies
W(z)*exp(W(z)) = z, and can thus be used to express solutions
of transcendental equations involving exponentials or logarithms.
The Lambert W function is also available in
Mathematica (as the ProductLog function), and in Maple.
}
\value{
Complex or real vector of solutions.
}
\references{Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), "On the Lambert
W Function", Advances in Computational Mathematics 5(4):329-359}
\author{Nici Schraudolph <schraudo at inf.ethz.ch> (original
version (c) 1998), Ben Bolker (R translation)
}
\note{
This implementation should return values within 2.5*eps of its
counterpart in Maple V, release 3 or later.  Please report any
discrepancies to the author or translator.
}
\examples{
curve(lambertW(x),from=0,to=10)
pvec = seq(-1,1,length=40)
m = outer(pvec,pvec,function(x,y)Re(lambertW(x+y*1i)))
persp(pvec,pvec,m,
num1 = uniroot(function(x) {x*exp(x)-1},lower=0,upper=1,tol=1e-9)
abs(lambertW(1)-num1\$root)<1e-9
}
\keyword{math}

--
620B Bartram Hall                            bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704