1 %% Copyright (C) 1998 by Nicol N. Schraudolph
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.
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
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/>.
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}.
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.
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.
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.
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>.
37 ## For further details, see:
39 ## Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), `On the Lambert
40 ## W Function', Advances in Computational Mathematics 5(4):329-359.
43 %% Author: Nicol N. Schraudolph <schraudo@inf.ethz.ch>
45 %% Created: 07 Aug 1998
46 %% Keywords: Lambert W Omega special transcendental function
48 function w = lambertw(b,z)
53 %% some error checking
57 if (any(round(real(b)) != b))
58 usage('branch number for lambertw must be integer')
63 %% series expansion about -1/e
65 % p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
70 % first-order version suffices:
72 w = (1 - 2*abs(b)).*sqrt(2*e*z + 2) - 1;
74 %% asymptotic expansion at 0 and Inf
76 v = log(z + ~(z | b)) + 2*pi*I*b;
79 %% choose strategy for initial guess
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;
92 t = f.*t./(p.*(w + f) - 0.5*(w + 2.0).*t./(w + f));
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))))
99 warning('iteration limit reached, result of lambertw may be inaccurate');