[R] Lambert's W function

[Stefan Böhringer]

On Tue, 2003-11-25 at 12:09, Robin Hankin wrote:
> Hello List
> 
> does anyone have an R function for the Lambert W function?  I need 
> complex arguments.
> 
> [the Lamert W function W(z) satisfies
> 
> W(z)*exp(W(z)) = z
> 
> but I could'nt even figure out how to use uniroot() for complex z]
> 
> 

There are several 'branches' of W(z) for complex arguments, since the
inverse of f(z) = z * exp(z) does not exist. Therefore you have to
constrain root-finding methods appropriately.

Here is an implementation for real valued arguments, if that is of any
help to you:

WE1 <- function(x) {
	if (1) {
	o <- optim(c(1),
		# optimize squared difference of z and log(t) + t
		function(t) { (x - (log(t[1]) + t[1]))**2 },
		# the gradient of the above function
		#function(t) { 2(log(t[1]) + t[1] - x)(1/t[1] - 1) },
		method="L-BFGS-B", lower=exp(-70),
			control = list(lmm=40, maxit=10000, factr=1e9));
	} else {
	o <- optim(c(1),
		# optimize squared difference of z and log(t) + t
		function(t) { (x - (log(t[1]) + t[1]))**2 },
		method="BFGS");
	}
	o$par
}

W2 <- function(x) {
	eps <- 10^-5;	# precision

	if (x == 0) {
		wnew <- 0
	} else {
		# initial guess
		wold <- if (x < 2.5) {
			(x + 4/3*x^2)/(1 + 7/3*x + 5/6*x^2)
		} else {
			log(x)
		};
		if (x < 0) print("error");
			lx <- log(x);
		wnew <- 2 * wold;
		if (x != Inf) {
			while (abs( (wnew - wold)/wold ) > eps) {
				wold <- wnew
				zn <- lx - wold - log(wold);
				y <- 2*(1 + wold)*((1 + wold) + 2/3*zn) - zn;
				wnew <- wold*(1 + (zn * y)/((1 + wold)*(y - zn)));
			}
		} else {
			wnew <- Inf;
		}
	}
	wnew
}
WE <- function(z) {
	if ( z < -1 ) { W2(exp(z)) }
	else { WE1(z) };
}
W <- function(z) { WE(log(z)) }


Stefan