Numerically Approximating the Principal Branch of Lambert's W Function

Numerically Approximating the Principal branch of the Complex Lambert's W function

We set forth now to try to approximate the infamous Lambert's W function, not only for real values but for complex values as well.

Let W denote the principal branch of Lambert's W function. Assume also that we will be working with values of theta in [-Pi,Pi], by defining our branch cut to be the negative real axis and our branch point to be 0.

Having seen a Maclaurin expansion of W in one of the previous articles, it is now interesting to examine the possibility of calculating the W for complex values outside its region of convergence: R = {z: |z| < 1/e}. Perhaps even calculate the non-principal branches of W as well.

The definition of the complex W is similar to the definition of W for the real case. Let us then try to solve for the complex number z that satisfies the equation:

z*e^z = w (1)

with complex w. Let:

w = R*e^i*Theta. (2) and
z = r*e^i*theta. (3)

Plugging (2) and (3) into (1), we get:

r*e^i*theta*e^{r*exp(i*theta)} = R*e^i*Theta.

The previous after some elementary calculations becomes:

r*e^r*cos(theta)*e^{i*(theta+r*sin(theta))} = R*e^i*Theta.

From the above, we immediately get a system of two equations:

R = r*e^r*cos(theta) (4)
Theta = theta + r*sin(theta) (5)

Our objective now is to effectively solve the above system for r and theta.

Note however, that we can immediately pick up a closed form expression for r from (4) via the very definition of W:

(4) => r = W(R*cos(theta))/cos(theta). (6)

Now we can use (6) along with (5) and reduce the system to the equation

(5)(6)=> Theta = theta + W(R*cos(theta))*tan(theta) (7).

Now, (7) can be solved using numerical methods, such as Newton's iteration.

The keen observer will remark here that equation (7) is circular: How can we determine solutions to equation (7) if, in fact we are trying to determine values for W?

The key point here is that we are assuming that the REAL equation x*e^x = a, CAN be solved numerically, via either W's maclaurin expansion if |x|<1/e, or via Newton's iteration if |x|>=1/e. And the instance of W we have above is essentially an instance of the REAL W. Or is it?

The reader here may want some justification as to the WHY Maple will be able to always solve it numerically. Let's investigate (7) some more:

First we write (7) in the form:

f(Theta,R,theta)= Theta - theta - W(R*cos(theta))*tan(theta) (8)

Main Claim:

(8) admits a real root in [0..Pi].

For this we need several preliminary Lemmas:

Lemma #1:

The REAL function g(x) = x*e^x is continuous and one-to one on the interval [0,a], for any a > 0.

Proof:

g is the product of two continuous and one-to-one functions in that interval.

Lemma #2:

The REAL W is continuous on the interval [0,a]

Proof:

This is theorem 3 on Spivak's page 220, with f = g from Lemma #1 above and f^-1 = W, since g(W(x)) = x, by the definition of the REAL W function.

Lemma #3

f of (8) has a removable discontinuity at theta = Pi/2.

Proof:

The bad term can be written as:

W(R*cos(theta))/cos(theta)*sin(theta)

We check the singularity around x=Pi/2:

lim_theta->Pi/2+W(R*cos(theta))/cos(theta), and
lim_theta->Pi/2-W(R*cos(theta))/cos(theta)

The above limits are essentially:

lim_x->0+W(R*x)/x, and
lim_x->0-W(R*x)/x

And these follow easily from the definition of Lambert's W function and the fact that W(0)=0:

W(R*x)*e^W(R*x)=R*x
=> W(R*x)/x = R/e^W(R*x).
=> lim_x->0+W(R*x)/x = R/e^W(R*0)= R. (9)

(Check that L'Hospital CANNOT be used here. :*)

Having established (9), we can redefine (8) as:

F(Theta,R,theta)= {Theta - theta - W(R*cos(theta))*tan(theta), theta is in [0,Pi/2) U (Pi/2,Pi]
Theta - Pi/2 - R, theta=Pi/2} (10)

and the result follows.

Now, we restrict the Real W so that F of (10) is always real-valued. This is achieved by forcing:
W(R*cos(theta)) >=-1, which then forces:
R*cos(theta) >=-1/e, which in turn forces
cos(theta) >= -1/[Re], which finally forces:
theta in [0..arccos(-1/[Re])].

Now let a = arccos(-1/[Re]). Note that a is always in (Pi/2, Pi). (11)

Lemma #4

F of (10) is continuous on [0..a].

Proof:

F is built using continuous functions on [0,Pi/2) U (Pi/2,a], so it is continuous there. Since by Lemma #3, it has a removable discontinuity at theta = Pi/2, F is continuous throughout [0,a].

We are now ready for the proof of the main claim: (8) admits a root in [0..a].

Proof:

F(Theta,R,0) =
=Theta-0-W(R*cos(0))*tan(0)
=Theta

If Theta = 0, we are done. theta=0 is a root.
If Theta > 0 then: (We are working with the upper half-plane only)

Let a = arccos(-1/[Re]).

F(Theta,R,a) =
=Theta-a-W(R*cos(a))*tan(a)
=Theta-a-W(-1/e)*tan(a)
=Theta-a+tan(a)

Now consider the function:
h(Theta,R) = Theta-a+tan(a)

h is continuous for fixed Theta in (Pi/2,Pi) as a function of R, is easily seen to be strictly decreasing as a function of R and additionally h(Theta,1/e) = Theta - Pi + tan(a) <= 0 (by 11). So h is negative for any Theta throughout (Pi/2,Pi). Therefore:

F(Theta,R,a) < 0.

Now we invoke continuity on F and the intermediate value theorem and this guarantees us a root in [0,a]. (!!)

We are done with the preliminaries. So let us summarize:

Assume we can solve the REAL equation: x*e^x = a, either via its Maclaurin expansion or via Newton's method. Call the real W: RW and our Complex W: CW.

Now then, given a complex number x+y*I:

if Im(w)=0 and w >= -1/e, then our case reduces to z=RW(w), so done.
if Im(w)<>0 then: {
if |z|<1/e then we approximate W via its Maclaurin expansion.
if |z|>=1/e then we proceed via equation (7), which in our case would be solving system (4) + (5) via numerical methods.
}

We begin with the code to solve the REAL case:
(A variant of the real case was provided by Robert Israel so we can have nested "fsolves". Note that it returns unevaluated unless the argument is numeric).

> restart;

> RW:=proc(a)
> local x;
> if type(a,numeric) then
> if a=evalf(-exp(-1)) then #special case W(-1/e)=-1
>   -1;
> elif a=evalf(exp(1)) then #special case W(e)=1
>   1;
> else
>   fsolve(x*exp(x)=a);
> fi;
> else
> 'procname(a)';
> fi;
> end:

Follows then some Maple code to approximate the Principal branch of the Complex Lambert's W. The observant reader will notice that several calculations are essentially redundant:

For example, since W(conjugate(z)) = conjugate(W(z)) we can restrict our calculations to the upper-half plane, so our Theta can range in [0..Pi].

We have to check the case z=-1/e manually, since Maple seems unable to solve the obvious equation: x*e^x = -1/e. (Which has x=-1 as a solution)

A detailed diagram on how the cases are checked is presented in the following picture:

The sky blue region which is the real case is checked first. Then the purple region which is the region of convergence is checked and finally the lower gray and upper purple gray half plane regions are checked. The lower half-plane falls back to the upper half-plane case, because of conjugation.

Here then, is the code for the Principal Branch of the Complex Lambert's W function:

> readlib(polar):

> CW:=proc(z)
> local R,Theta,r,theta,rez,imz;
> rez:=evalf(Re(z));imz:=evalf(Im(z));
> if (imz=0) and (evalf(-exp(-1))<=rez) then
> RW(z); #REAL case!!
> else
> if evalf(abs(z))<evalf(exp(-1)) then #inside region of convergence.
>   sum((-1)^(n-1)*n^(n-1)*z^n/n!,n=1..30); See expansion of W
> else #outside region of convergence
>   if (imz<0) then #call recursively for conjugates.
>    conjugate(CW(conjugate(z))); #See Lemma #3 on Infinite Exponentials
>   else #all other cases, see analysis above!!
>    R:=evalf(op(1,polar(z)));   #Turn z into polar form
>    Theta:=evalf(op(2,polar(z)));
>    theta:=fsolve(Theta=theta+RW(R*cos(theta))*tan(theta),theta,0..evalf(Pi));
>    if theta=evalf(Pi/2) then
>     r:=R;
>    else
>     r:=RW(R*cos(theta))/cos(theta);
>    fi;
>    evalf(r*cos(theta)+r*I*sin(theta));
>   fi;
> fi;
> fi;
> end:

Finally we compare our procedure to the W function built in Maple:

>W:=LambertW;

> evalf(W(3+4*I));
1.281561806 + .5330952220 I

> evalf(CW(evalf(3+4*I)));
1.281561806 + .5330952220 I

> evalf(W(-exp(-1)+0.001+0.1*I));
-.5031781151 + .3773368931 I

> evalf(CW(evalf(-exp(-1)+0.001+0.1*I)));
-.5031781147 + .3773368930 I

evalf(W(-20));
1.908615873 + 2.269938354 I

> evalf(CW(-20));
1.908615873 + 2.269938354 I

evalf(W(0.02+0.01*I));
.01970316808 + .009615880874 I

> evalf(CW(evalf(0.02+0.01*I)));
.01970316808 + .009615880874 I

Obviously the code for CW above can be optimized further to avoid complexes altogether, to speed up certain calculations and to port this to Pascal, C or FORTRAN.

Note also how we take care of the special case theta=Pi/2, by modifying the r function based on the limit R of equation (10) to avoid stepping onto Pi/2 into the code for CW.

You'd have to take care of all the special real values by using conventional real comparisons, which Maple is able to do here in very simple ways: For example, in most programming languages one never checks whether r=some real value directly, since this is guaranteed to be almost always false. Instead, if a is a known real value and x is an unknown real value, one always checks equality via: |x-a| < epsilon, where epsilon is a given small number.

The specific implementation details are left as an exercise for the reader. An alternative algorithm for calculating complex branches of W, appears here.

Back to Mathematics