How to solve Lambert W function?

Specific Domains Numerics
1

I have this equation

Nₜ₊₁ = Nₜ exp[μ(1-Nₜ/K)]

I would like to solve it the other way round: knowing Nₜ₊₁, I would like to estimate Nₜ. This is in the form
y = x exp(a(1−x/k))
From this post looks like the solution is not so straightforward since it requires to solve the Lambert function W, which involves complex numbers.

Is there a way to solve this function with julia?

Thank you

1 Like
2

GSL.jl provides an interface to the GNU Scientific Library, which has Lambert’s W.

2 Likes
3

I suggest to use a root finder, eg Roots.jl or IntervalRootFinding.jl.

2 Likes
4

Newton’s method for f(w) = we^w - x = 0 would certainly be easy to implement here. Even higher-order iterations for the Lambert W function were given by Fritsch et al (1973) (though you should be careful to look at the math only, and not the source code, to avoid the restrictive ACM copyright terms), as well as an initial guess w \approx \frac{x + \frac{4}{3} x^2}{1 + \frac{7}{3}x + \frac{5}{6}x^2} based on a Padé approximant for |x| \lesssim 0.7385 and w \approx \log(x) - 24\frac{\log^2(x) + 2\log(x) - 3}{7\log^2(x) + 58\log(x) +127} for larger |x|.

7 Likes
5

Is there something wrong with the LambertW.jl package? LambertW | Julia Observer

2 Likes
6

Nothing AFAIK — I just forgot about it. (It might be possible to speed it up using Fritsch’s algorithm, as LambertW.jl currently seems to use the slower Halley’s iteration.)

5 Likes
7

Maybe that can be integrated into the SpecialFunctions package?

4 Likes
8

I tend to agree: it is a little peculiar that there is a separate package for the Lambert W function.

Anyway, I was looking into pipe friction some years ago, and found a paper using the Lambert W function to solve the model. The fact that Julia had the LambertW.jl package was actually what made me switch from Python to Julia…

5 Likes
9

The Lambert W function does not necessarily involve complex numbers.

1 Like
10

LambertW.jl could certainly be speeded up. You’d have to benchmark using a more accurate initial condition, it might be more costly than a few more iterations of the root finder.

There was talk of integrating it into SpecialFunctions.jl a few years ago. I don’t recall why it stalled. Julia is not node.js, but it’s not fortran either. It’s not obvious to me that special functions need to be collected in one place. But, I’m not against it, either.

There is an implementation of the Lambert W function in Python’s mpmath, btw.

1 Like
11

I’m looking at it and will see if it can be sped up at all. The good news is it already should be a bit faster in 1.6 since exp has gotten faster and that looks like it should be where most of the time is spent.

1 Like
12

(The Fritsch recurrence uses log, not exp.)

1 Like
13

Sorry for being unclear. The halley iteration that it currently uses is exp based, so that is what’s getting a speedup in 1.6.

1 Like
14

As a quick test, the following lw(x) function (which only supports nonnegative real arguments) is about 3x faster in Float64 precision than LambertW.lambertw(x) on my machine with Julia 1.5:

lw(x) = _lw(float(x))
function _lw(x::T) where {T<:Real}
    x < 0 && error("only x ≥ 0 supported")
    
    # Based on the algorithm described in 
    #    F. N. Fritsch, R. E. Shafer, and W. P. Crowley, “Solution of the transcendental equation we w = x,”
    #    Commun. ACM, vol. 16, no. 2, pp. 123–124, Feb. 1973, doi: 10.1145/361952.361970. 
    # (without looking at the accompanying source code):
    logx = log(x)
    if x < 0.7385
        w = x * @evalpoly(x, 6, 8) / @evalpoly(x, 6, 14, 5)
    else
        w = logx + @evalpoly(logx, 72, -48, -24) / @evalpoly(logx, 127, 58, 7)
    end
    ε = sqrt(eps(typeof(w))) # quite conservative, since we should have quartic convergence
    while true
        z = logx - log(w) - w
        t1 = 1+w
        t2 = 2t1*(t1+z*(T(2)/T(3))) - z
        e = z * t2 / (t1 * (t2 - z))
        w += e*w
        abs(e) < ε && break
    end
    return w
end

and it could be made faster by sharpening the convergence test — it seems like it does one more iteration than it needs in many cases. (In Float64 precision, it invariably seems to require only one iteration to reach machine precision, so it is possible that the loop could be omitted entirely in that case.)

11 Likes
15

Thank you all!
The problem is I am not a mathematician, thus I find all this a bit overwhelming. In Particular, I don’t understand how do I get x from _lw (from where the parameters come from? z = logx - log(w) - w does not look y = x exp(a(1−x/k)) to me…).
I looked at lambertW.jl, but it looks like this function returns a number ( lambertw(10)~1.7). How do I use it to get Nₜ?

16

To solve this for x, rewrite it as -ae^{-a} y/k = (-ax/k) e^{-ax/k}. Then -ax/k = \mathrm{W}(-ae^{-a} y/k) where \mathrm{W} is the Lambert W function, so x = -(k/a)\mathrm{W}(-ae^{-a} y/k).

In Julia with the LambertW.jl package, for example, x = -(k/a)*lambertw(-a*exp(-a)*y/k).

4 Likes
17

you have to check that you are on the right branch. I think -ax/k>-1/e

1 Like
18

I ran the following:

function ricky(x, k, r, I)
  for i in 1:I
    q = 1-x/k
    y = x * exp(r*q)
    x = y
  end
  return y
end
y = ricky(50000, 2.2*10^7, 0.47, 2)
# a given population increases from 50 000 to 79 914.30 between to and t1
# now, let's invert it
using LambertW
x = -(k/r)*lambertw(-r*exp(-r)*y/k)
julia> x
49999.99999999999

which looks pretty good to me…
Thank you!

3 Likes