lambertw_yp

 LUKE - LAMBERTW Lambert's W function.

    LAMBERTW Lambert's W function. 
   W = LAMBERTW(X) is the solution to w*exp(w) = x.
   W = LAMBERTW(K,X) is the K-th branch of this multi-valued function. Z may be a 
   complex scalar or array.  For real Z, the result is real on
   the principal branch for Z >= -1/e.

   W = 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.

    Input:

        - b: data matrix [n,m]
        - z: data matrix [n,m]
       - mode: (0) full numerical approach (fast, any matlab version), (1) using Symbolic Math Toolbox is exist and asked (quite slow)
       
          (default = '0')

    Output:
    
        - w: Lambert's W function [n,m]

by Y.PEYSSON CEA/IRFM <yves.peysson@cea.fr>

0001 function w = lambertw_yp(b,z,mode)
0002 % LUKE - LAMBERTW Lambert's W function.
0003 %
0004 %    LAMBERTW Lambert's W function.
0005 %   W = LAMBERTW(X) is the solution to w*exp(w) = x.
0006 %   W = LAMBERTW(K,X) is the K-th branch of this multi-valued function. Z may be a
0007 %   complex scalar or array.  For real Z, the result is real on
0008 %   the principal branch for Z >= -1/e.
0009 %
0010 %   W = LAMBERTW(B,Z) specifies which branch of the Lambert
0011 %   W-Function to compute.  If Z is an array, B may be either an
0012 %   integer array of the same size as Z or an integer scalar.
0013 %   If Z is a scalar, B may be an array of any size.
0014 %
0015 %   The algorithm uses series approximations as initializations
0016 %   and Halley's method as developed in Corless, Gonnet, Hare,
0017 %   Jeffrey, Knuth, "On the Lambert W Function", Advances in
0018 %   Computational Mathematics, volume 5, 1996, pp. 329-359.
0019 %
0020 %    Input:
0021 %
0022 %        - b: data matrix [n,m]
0023 %        - z: data matrix [n,m]
0024 %       - mode: (0) full numerical approach (fast, any matlab version), (1) using Symbolic Math Toolbox is exist and asked (quite slow)
0025 %
0026 %          (default = '0')
0027 %
0028 %    Output:
0029 %
0030 %        - w: Lambert's W function [n,m]
0031 %
0032 %by Y.PEYSSON CEA/IRFM <yves.peysson@cea.fr>
0033 %
0034 [dummy,dummy,flag_toolbox] = matversion_yp('Symbolic Math Toolbox');
0035 %
0036 if flag_toolbox == 0 || nargin < 3, mode = 0;end
0037 %
0038 if nargin == 1,
0039     z = b;
0040     b = 0;
0041 end
0042 %
0043 if mode == 0,
0044     %
0045     % Use asymptotic expansion w = log(z) - log(log(z)) for most z
0046     %
0047     tmp = log(z + (z == 0)) + i*b*6.28318530717958648;
0048     w = tmp - log(tmp + (tmp == 0));
0049     %
0050     % For b = 0, use a series expansion when close to the branch point
0051     %
0052     k = find(b == 0 & abs(z + 0.3678794411714423216) <= 1.5);
0053     tmp = sqrt(5.43656365691809047*z + 2) - 1 + i*b*6.28318530717958648;
0054     w(k) = tmp(k);
0055     %
0056     for k = 1:36
0057         %
0058         % Converge with Halley's iterations, about 5 iterations satisfies the tolerance for most z
0059         %
0060         c1 = exp(w);
0061         c2 = w.*c1 - z;
0062         w1 = w + (w ~= -1);
0063         dw = c2./(c1.*w1 - ((w + 2).*c2./(2*w1)));
0064         w = w - dw;
0065         %
0066         if all(abs(dw) < 0.7e-16*(2+abs(w)))
0067             break;
0068         end
0069     end
0070     %
0071     % Specially handle z = 0
0072     %
0073     w(b ~= 0 & z == 0) = -inf;
0074     %
0075 else
0076     w = lambertw(b,z); 
0077 end
0078