octave_packages/specfun-1.1.0/lambertw.m

   1 %% Copyright (C) 1998 by Nicol N. Schraudolph
   2 %%
   3 %% This program is free software; you can redistribute and/or
   4 %% modify it under the terms of the GNU General Public
   5 %% License as published by the Free Software Foundation;
   6 %% either version 3, or (at your option) any later version.
   7 %%
   8 %% This program is distributed in the hope that it will be
   9 %% useful, but WITHOUT ANY WARRANTY; without even the implied
  10 %% warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11 %% PURPOSE.  See the GNU General Public License for more
  12 %% details.
  13 %%
  14 %% You should have received a copy of the GNU General Public
  15 %% License along with this software; see the file COPYING.  If not,
  16 %% see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {@var{x} = } lambertw (@var{z})
  20 ## @deftypefnx {Function File} {@var{x} = } lambertw (@var{z}, @var{n})
  21 ## Compute the Lambert W function of @var{z}.
  22 ##
  23 ## This function satisfies W(z).*exp(W(z)) = z, and can thus be used to express
  24 ## solutions of transcendental equations involving exponentials or logarithms.
  25 ##
  26 ## @var{n} must be integer, and specifies the branch of W to be computed;
  27 ## W(z) is a shorthand for W(0,z), the principal branch.  Branches
  28 ## 0 and -1 are the only ones that can take on non-complex values.
  29 ##
  30 ## If either @var{n} or @var{z} are non-scalar, the function is mapped to each
  31 ## element; both may be non-scalar provided their dimensions agree.
  32 ##
  33 ## This implementation should return values within 2.5*eps of its
  34 ## counterpart in Maple V, release 3 or later.  Please report any
  35 ## discrepancies to the author, Nici Schraudolph <schraudo@@inf.ethz.ch>.
  36 ##
  37 ## For further details, see:
  38 ##
  39 ## Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), `On the Lambert
  40 ## W Function', Advances in Computational Mathematics 5(4):329-359.
  41 ## @end deftypefn
  42
  43 %% Author:   Nicol N. Schraudolph <schraudo@inf.ethz.ch>
  44 %% Version:  1.0
  45 %% Created:  07 Aug 1998
  46 %% Keywords: Lambert W Omega special transcendental function
  47
  48 function w = lambertw(b,z)
  49     if (nargin == 1)
  50         z = b;
  51         b = 0;
  52     else
  53         %% some error checking
  54         if (nargin != 2)
  55             print_usage;
  56         else
  57             if (any(round(real(b)) != b))
  58                 usage('branch number for lambertw must be integer')
  59             end
  60         end
  61     end
  62
  63     %% series expansion about -1/e
  64     %
  65     % p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
  66     % w = (11/72)*p;
  67     % w = (w - 1/3).*p;
  68     % w = (w + 1).*p - 1
  69     %
  70     % first-order version suffices:
  71     %
  72     w = (1 - 2*abs(b)).*sqrt(2*e*z + 2) - 1;
  73
  74     %% asymptotic expansion at 0 and Inf
  75     %
  76     v = log(z + ~(z | b)) + 2*pi*I*b;
  77     v = v - log(v + ~v);
  78
  79     %% choose strategy for initial guess
  80     %
  81     c = abs(z + 1/e);
  82     c = (c > 1.45 - 1.1*abs(b));
  83     c = c | (b.*imag(z) > 0) | (~imag(z) & (b == 1));
  84     w = (1 - c).*w + c.*v;
  85
  86     %% Halley iteration
  87     %
  88     for n = 1:10
  89         p = exp(w);
  90         t = w.*p - z;
  91         f = (w != -1);
  92         t = f.*t./(p.*(w + f) - 0.5*(w + 2.0).*t./(w + f));
  93         w = w - t;
  94         if (abs(real(t)) < (2.48*eps)*(1.0 + abs(real(w)))
  95             && abs(imag(t)) < (2.48*eps)*(1.0 + abs(imag(w))))
  96             return
  97         end
  98     end
  99     warning('iteration limit reached, result of lambertw may be inaccurate');
 100 end
 101