In the previous post we wanted to find a value of *n* such that *f*(*n*) = 10^{12} where

and we took a rough guess *n* = 200. Turns out *f*(200) ≈ 4 × 10^{12} and that was good enough for our purposes. But what if we wanted to solve f(*n*) = *x* accurately?

We will work our way up to solving *f*(*n*) = 10^{12} exactly using the Lambert *W* function.

## Lambert’s *W* function

If you can solve one equation involving exponential functions you can bootstrap your solution solve a lot more equations.

The Lambert *W* function is defined to be the function *W*(*x*) that maps each *x* to a solution of the equation

*w* exp(*w*) = *x*.

This function is implemented Python under `scipy.special.lambertw`

. Let’s see how we could use it solve similar equations to the one that define it.

## Basic bootstrapping

Given a constant *c* we can solve

*cw* exp(w) = *x*

simply by solving

*w* exp(*w*) = *x*/*c*,

which says w* =* *W*(*x*/*c*).

And we can solve

*w* exp(*cw*) = *x*

by setting *y* = *cw* and solving

*y* exp(*y*) = *cx.*

and so *cw* = *W*(*cx*) and *w* = *W*(*cx*)/*c*.

Combining the two approaches above shows that the solution to

*aw* exp(*bw*) = *x*

is

*w* = *W*(*bx*/*a*) / *b*.

We can solve

exp(*w*) / *w* = *x*

by taking the reciprocal of both sides and applying the result above with *a* = 1 and *b* = -1.

More generally, we can solve

*w*^{c} exp(w) = *x*

by raising both sides to the power 1/*c*; the reciprocal the special case *c* = 1.

Similarly, we can solve

*w* exp(*w*^{c}) = *x*

by setting *y* = *w*^{c} , changing the problem to

*y*^{−1/c} exp(*y*) = *x*.

## Putting it all together

We’ve now found how to deal with constants and exponents on *w*, inside and outside the exponential function. We now have all the elements to solve our original problem.

We can solve the general equation

with

The equation at the top of the post corresponds to *a* = 1/√48, *b* = -1, *c* = π √(2/3), and *d* = 1/2.

We can code this up in Python as follows.

from scipy.special import lambertw def soln(a, b, c, d, x): t = (x/a)**(d/b) *c*d/b return (lambertw(t)*b/(c*d))**(1/d)

We can verify that our solution really is a solution by running it though

from numpy import exp, pi def f(a, b, c, d, w): return a*w**b * exp(c*w**d)

to make sure we get our *x* back.

When we run

a = 0.25*3**-0.5 b = -1 c = pi*(2/3)**0.5 d = 0.5 x = 10**12 w = soln(a, b, c, d, x) print(f(a, b, c, d, w))

we do indeed get *x* back.

a = 0.25*3**-0.5 b = -1 c = pi*(2/3)**0.5 d = 0.5 w = soln(a, b, c, d, 10) print(f(a, b, c, d, w))

When we take a look at the solution *w*, we get 1.443377079584187e-13. In other words, we get a number near zero. But our initial guess was *w* = 200. What went wrong?

Nothing went wrong. We got a solution to our equation. It just wasn’t the solution we expected.

The Lambert *W* function has multiple branches when viewed as a function on complex numbers. It turns out the solution we were expecting is on the branch SciPy denotes with `k = -1`

. If we add this to the call to `lambertw`

above we get the solution 183.85249773880142 which is more in line with what we expected.