# How to solve Lambert W function?

**URL:** <https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538>\
**Category:** Numerics\
**Tags:** question\
**Created:** [December 28, 2020, 9:04pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538 "2020-12-28T21:04:03Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 28, 2020, 9:04pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/1 "2020-12-28T21:04:03Z")

</div>

I have this equation

```julia
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](https://math.stackexchange.com/questions/3964810/how-to-solve-equation-for-a-variable-that-is-inside-and-outside-an-exponent) 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

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [December 28, 2020, 9:24pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/2 "2020-12-28T21:24:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [December 28, 2020, 10:25pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/3 "2020-12-28T22:25:30Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 28, 2020, 10:41pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/4 "2020-12-28T22:41:09Z")

</div>

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)](https://doi.org/10.1145/361952.361970) (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|.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [December 28, 2020, 10:49pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/5 "2020-12-28T22:49:43Z")

</div>

Is there something wrong with the `LambertW.jl` package? [LambertW | Julia Observer](https://juliaobserver.com/packages/LambertW)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 28, 2020, 10:54pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/6 "2020-12-28T22:54:05Z")

</div>

> [@BLI](#):
>
> Is there something wrong with the `LambertW.jl` package?

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.)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [December 28, 2020, 10:56pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/7 "2020-12-28T22:56:31Z")

</div>

Maybe that can be integrated into the SpecialFunctions package?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [December 28, 2020, 11:00pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/8 "2020-12-28T23:00:58Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [December 28, 2020, 11:36pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/9 "2020-12-28T23:36:06Z")

</div>

> [@Luigi\_Marongiu](#):
>
> since it requires to solve the Lambert function W, which involves complex numbers.

The Lambert W function does not necessarily involve complex numbers.

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [December 28, 2020, 11:48pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/10 "2020-12-28T23:48:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 29, 2020, 12:05am UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/11 "2020-12-29T00:05:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 29, 2020, 1:14am UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/12 "2020-12-29T01:14:13Z")

</div>

> [@Oscar\_Smith](#):
>
> 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.

(The Fritsch recurrence uses `log`, not `exp`.)

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 29, 2020, 1:16am UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/13 "2020-12-29T01:16:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 29, 2020, 1:37am UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/14 "2020-12-29T01:37:08Z")

</div>

> [@jlapeyre](#):
>
> 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.

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:

```julia
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.)

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 29, 2020, 9:52am UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/15 "2020-12-29T09:52:53Z")

</div>

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.](https://www.juliaobserver.com/packages/LambertW)jl, but it looks like this function returns a number ( `lambertw(10)~1.7`). How do I use it to get Nₜ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 29, 2020, 2:03pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/16 "2020-12-29T14:03:47Z")

</div>

> [@Luigi\_Marongiu](#):
>
> 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.

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)`.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 29, 2020, 3:34pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/17 "2020-12-29T15:34:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 29, 2020, 3:48pm UTC](https://discourse.julialang.org/t/how-to-solve-lambert-w-function/52538/18 "2020-12-29T15:48:07Z")

</div>

I ran the following:

```julia
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!