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
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
