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 θ in [-π,π], 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*ez = w (1)
with complex w. Let:
w = R*ei*Θ. (2) and
z = r*ei*θ. (3)
Plugging (2) and (3) into (1), we get:
r*ei*θ*er*exp(i*θ) = R*ei*Θ.
The previous after some elementary calculations becomes:
r*er*cos(θ)*ei*(θ+r*sin(θ)) = R*ei*Θ.
From the above, we immediately get a system of two equations:
R = r*er*cos(θ) (4)
Θ = θ + r*sin(θ) (5)
Our objective now is to effectively solve the above system for r and θ.
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(θ))/cos(θ). (6)
Now we can use (6) along with (5) and reduce the system to the equation
(5)(6)=> Θ = θ + W(R*cos(θ))*tan(θ) (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*ex = a, can be solved numerically, via either W's 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(Θ,R,θ)= Θ - θ -W(R*cos(θ))*tan(θ) (8)
Main Claim:
(8) admits a real root in [0..π].
For this we need several preliminary Lemmas:
Lemma #1:
The real function g(x) = x*ex 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 θ = π/2.
Proof:
The bad term can be written as:
W(R*cos(θ))/cos(θ)*sin(θ)
We check the singularity around x=π/2:
limθ->π/2+W(R*cos(θ))/cos(θ), and limθ->π/2-W(R*cos(θ))/cos(θ)
The above limits are essentially:
limx->0+W(R*x)/x, and limx->0-W(R*x)/x
These follow easily from the definition of Lambert's W function and the fact that W(0)=0:
W(R*x)*eW(R*x)=R*x
=> W(R*x)/x = R/eW(R*x).
=> limx->0+W(R*x)/x = R/eW(R*0)= R. (9)
(Check that L'Hospital cannot be used here.
Having established (9), we can redefine (8) as:
F(Θ,R,θ)= {Θ - θ - W(R*cos(θ))*tan(θ), θ is in [0,π/2) U (π/2,π], Θ - π/2 - R, θ=π/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(θ)) >=-1, which then forces:
R*cos(θ) >=-1/e, which in turn forces
cos(θ) >= -1/[Re], which finally forces:
θ in [0..arccos(-1/[Re])].
Now let a=arccos(-1/[Re]). Note that a is always in (π/2, π). (11)
Lemma #4
F of (10) is continuous in [0..a].
Proof:
F is built using continuous functions on [0,π/2) U (π/2,a], so it is continuous there. Since by Lemma #3, it has a removable discontinuity at θ = π/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(Θ,R,0) =
=Θ-0-W(R*cos(0))*tan(0)=Θ
If Θ = 0, we are done. θ=0 is a root. If Θ > 0 then (we are working with the upper half-plane only) let a = arccos(-1/[Re]).
F(Θ,R,a) =
=Θ-a-W(R*cos(a))*tan(a)
=Θ-a-W(-1/e)*tan(a)
=Θ-a+tan(a)
Now consider the function:
h(Θ,R) = Θ-a+tan(a)
h is continuous for fixed Θ in (π/2,π) as a function of R, is easily seen to be strictly decreasing as a function of R and additionally h(Θ,1/e) = Θ - π + tan(a) ≤0 (by 11). So h is negative for any Θ throughout (π/2,π). Therefore F(Θ,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 preliminaries. Let us summarize. Assume we can solve the real equation: x*ex = a, either via its expansion or via Newton's method. Call the real W, RW and the Complex W, CW.
Given a complex number x+y*i, if Im(w)=0 and w >=-1/e, then our case reduces to z=RW(w), so we are 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 some Maple code to approximate the principal branch of the Complex Lambert's W. The 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 Θ can range in [0..π]. We have to check the case z=-1/e manually, since Maple seems unable to solve the obvious equation: x*ex = -1/e, which has x=-1 as a solution.
A detailed diagram on how the cases are checked is presented in the following figure.
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:
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
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 how we take care of the special case θ=π/2, by modifying the r function based on the limit R of equation (10) to avoid stepping onto π/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| < ε, where ε is a given small positive number.
The specific implementation details are left as an exercise for the reader[1].