Solving equations with Lambert’s W function
In the previous post we wanted to find a value of n such that f(n) = 1012 where
and we took a rough guess n = 200. Turns out f(200) 4 * 1012 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) = 1012 exactly using the Lambert W function.
Lambert's W functionIf 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 bootstrappingGiven 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
wc 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(wc) = x
by setting y = wc , changing the problem to
y-1/c exp(y) = x.
Putting it all togehterWe'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.
More Lambert W postsThe post Solving equations with Lambert's W function first appeared on John D. Cook.