Name

specfun_lambertw — The Lambert W-Function.

Calling Sequence

   w = specfun_lambertw ( z )
   w = specfun_lambertw ( b , z )
   
   

Parameters

b :

a matrix of doubles (default b = 0)

z :

a matrix of doubles

w :

a matrix of doubles

Description

W = specfun_lambertw(Z) computes the principal value of the Lambert W-Function, the solution of Z = W*exp(W). Z may be a complex scalar or array. For real Z, the result is real on the principal branch for Z >= -1/e.

W = specfun_lambertw(B,Z) specifies which branch of the Lambert W-Function to compute. If Z is an array, B may be either an integer array of the same size as Z or an integer scalar. If Z is a scalar, B may be an array of any size.

The algorithm uses series approximations as initializations and Halley's method as developed in Corless, Gonnet, Hare, Jeffrey, Knuth, "On the Lambert W Function", Advances in Computational Mathematics, volume 5, 1996, pp. 329-359.

Examples

c = specfun_lambertw([0 -exp(-1); %pi 1])
e = [
0                     -1.
1.07365819479614921    0.56714329040978384
]

// The following is an Inversion Error Test.
// Ideally, lambertw(b,z)*exp(lambertw(b,z)) = z for any complex z
// and any integer branch index b, but this is limited by machine
// precision.  The inversion error |lambertw(b,z)*exp(lambertw(b,z)) - z|
// is small but worth minding.

// Experimentation finds that the error is usually on the order of
// |z|*1e-16 on the principal branch.  This test computes the inversion
// error over the square [-10,10]x[-10,10] in the complex plane, large
// enough to characterize the error away from the branch points at
// z = 0 and -1/e.

// Use NxN points to sample the complex plane
N = 81;
// Sample in the square [-R,R]x[-R,R]
R = 10;
x = linspace(-R,R,N);
y = linspace(-R,R,N);
[xx,yy] = meshgrid(x,y);
z = xx + %i*yy;
i = 1;
j = 1;
msubp = 3;
nsubp = 3;
for b = -4:4
w = specfun_lambertw(b,z);
InvError = abs(w.*exp(w) - z);
mprintf("Largest error for b = %2d:  %.2e\n",b,max(InvError));
p = (i-1)*nsubp + j;
subplot(msubp,nsubp,p);
h = gcf();
cmap = graycolormap (10);
h. color_map = cmap ;
surf(x,y,InvError');
xtitle("Inversion error for W_"+string(b)+"(z)" , "Re z" , "Im z" , "Error" );
j = j + 1;
if ( j > nsubp ) then
j = 1;
i = i + 1;
end
end

   

Authors

Pascal Getreuer 2005-2006
Modified by Didier Clamond, 2005
Michael Baudin, DIGITEO, 2010 (Scilab port)

Bibliography

http://www.mathworks.com/matlabcentral/fileexchange/6909