<< specfun_gammainc specfun specfun_log1p >>

specfun >> specfun > specfun_lambertw

specfun_lambertw

The Lambert W-Function.

Calling Sequence

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

Parameters

z :

a matrix of doubles

b :

a matrix of floating point integers, the index of the branch of the multi-valued function (default b = 0).

w :

a matrix of doubles

Description

If b=0, this function computes the principal value of the Lambert W-Function, the solution of the equation

z may be a complex scalar or array. For real z, the result is real on the principal branch for z >= -1/e.

The parameter b 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.

Any optional argument equal to the empty matrix [] is replaced by its default value.

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.

The calling sequence specfun_lambertw ( z , b ) is different from Matlab/lambertw: the two arguments z and b are switched. The current choice is kept in Scilab, because this keeps the management of the optionnal input argument b consistent.

Examples

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

// Compute Wk(1), for k=-4,-3,...,4
w = specfun_lambertw(1,(-4:4)')
e = [
-3.162952738804083896D+00-%i*2.342774750375521364D+01;
-2.853581755409037690D+00-%i*1.711353553941214756D+01;
-2.401585104868003029D+00-%i*1.077629951611507053D+01;
-1.533913319793574370D+00-%i*4.375185153061898369D+00;
5.671432904097838401D-01;
-1.533913319793574370D+00+%i*4.375185153061898369D+00;
-2.401585104868003029D+00+%i*1.077629951611507053D+01;
-2.853581755409037690D+00+%i*1.711353553941214756D+01;
-3.162952738804083896D+00+%i*2.342774750375521364D+01
];

// Produce a plot of |W0(z)|.
x = linspace(-6,6,51);
y = linspace(-6,6,51);
[x,y] = meshgrid(x,y);
w = specfun_lambertw(x + %i*y);
surf(x,y,abs(w));
xlabel("Re z");
ylabel("Im z");
title("|W_0(z)|");

// The following is an Inversion Error Test.
// Ideally, lambertw(z,b)*exp(lambertw(z,b)) = z for any complex z
// and any integer branch index b, but this is limited by machine
// precision.  The inversion error |lambertw(z,b)*exp(lambertw(z,b)) - 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(z,b);
    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

Bibliography

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

Lambert W Function, Pascal Getreuer, http://www.math.ucla.edu/~getreuer/lambertw.html

Corless, Gonnet, Hare, Jeffrey, Knuth, "On the Lambert W Function", Advances in Computational Mathematics, volume 5, 1996, pp. 329-359.


Report an issue
<< specfun_gammainc specfun specfun_log1p >>