specfun_lambertw — The Lambert W-Function.
w = specfun_lambertw ( z ) w = specfun_lambertw ( b , z )
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.
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