TR-00-20.mws

• Introduction
  • History
  • Lambert's Trinomial Equation (in Euler's format)
• Basic Facts about W
  • Applications
    • Iteration of the exponential function (Euler, Carlsson, Hayes, Wright, ...)
  • Complex Analysis
    • Branches
      • The Quadratrix of Hippias
      • An Animation of a few values of W(k,z) for varying z and k (recondite, says David)
    • Riemann Surface
      • The Animated Riemann Surface for The Lambert W function
      • A proof that the computed surface is correct
  • Calculus
    • Derivatives
    • Series
    • Integration
  • Computation
• A surprising integral
• The Wright
  omega; function
  • Applications
  • Complex Analysis
    • Branches
    • Riemann Surface
      • Actual plot of the Riemann Surface for
        omega;.
      • Proof that this computation works (incomplete)
  • Calculus
    • Derivatives
    • Series
    • Integration
  • Computation
• Concluding Remarks

---

Still more fun results on
the Lambert W function

Robert M. Corless and David J. Jeffrey
Department of Applied Mathematics
University of Western Ontario

Introduction

This worksheet explores some recent results related to the function W(x), which satisfies

> restart;

> W(x)*exp(W(x)) = x;

In fact, Maple knows this function rather well, and names it LambertW, following the paper ``On the Lambert W function'', by Corless, Gonnet, Hare, Jeffrey, and Knuth (Adv. Comp. Math. 1996); to save typing its real name all the time, we use

> alias( W = LambertW );

This means now that any instance of W that occurs in this worksheet will use the short notation W instead of the long notation LambertW.

> solve( y*exp(y) - z, y );

> _EnvAllSolutions := true;

> solve( y*exp(y) - z, y );

>

The usual reaction to this answer is another question: what on earth is W(x)?

In this worksheet we will see some parts of the answer, and some of the history and applications of this function. The goal of this exposition is that, at the end, you feel comfortable with W(x) and are happy with it as an answer. Of course, to get to that point you will have to do some work, but luckily it is all rather pleasant.

History

The history of the function goes back to J. H. Lambert (1728-1777). This worksheet is not the right medium to discuss the life of Lambert, but a short note is appropriate.

Lambert was born the son of a tailor, and was expected by his father to continue in that profession. His early fight for his education is a remarkable story---he had to sell drawings and writings to his classmates to buy candles for night study, for example---but eventually his talents were recognized and he got a position as tutor in a house which had a decent library. He then was able to educate himself, and went on to produce fundamental discoveries in cartography (the Lambert projection is still in use), hygrometry, pyrometry, statistics, philosophy (where he is actually more famous than as a mathematician), and pure mathematics. He is most noted as being the first person to prove that is irrational, which was an important step in proving that the classical problem of squaring the circle was impossible by straightedge and compass.

Lambert's Trinomial Equation (in Euler's format)

Consider the `trinomial equation' (using Euler's formulation, not Lambert's)

> restart;

> Euler_trinomial := x^alpha - x^beta = (alpha-beta)*v*x^(alpha+beta);

> isolate( Euler_trinomial, v);

> simplify(%,symbolic);

> map( series, %, x=1, 6):

> map(factor,%);

Now reverse that series, to get a series for x in terms of v.

> solve(%, x):

> map(factor, %);

>

Lambert found that series as a solution to the trinomial equation, by hand calculation of the first few terms (his proof of the general case was, I believe, lost). His calculation was essentially similar to the above, except he did not carry so many terms, and used a more clever approach than the general series reversion techniques used by `solve' above; indeed his techniques proved a forerunner to the Lagrange Inversion Theorem which is now the classical way for reversion of series.

Lambert at about the same time discovered that the series for x raised to an arbitrary power n had a very similar form:

> x := %:

> series( x^n, v, 5):

> map(factor,%);

>

Lambert thought the pattern beautiful, and so simple that it must have a finite formula as a closed form. (I don't know of any, though, and I rather expect that it is in fact impossible).

Euler came somewhat closer to proving that this simple pattern for the coefficients continues for all terms of the series; however, our translation of his Latin paper is as yet incomplete (Matt Davison has done most of it so far, and a student will complete the process this summer). Euler at least showed that the function

> phi(alpha, beta, n, v);

defined by the infinite series, satisfies

> phi(alpha, beta, -alpha, v) - phi(alpha, beta, -beta, v) = (alpha-beta)*v;

This does not mean that

> phi(alpha,beta,n,v) = phi(alpha,beta,1,v)^n;

>

and this gap remained (unnoticed? if indeed Euler didn't close it himself in some part of the paper as yet untranslated). I do not know who first established that this series really does continue in the same way, but at the very latest, G. Raney in 1959 provided a combinatorial proof equivalent to the above identity (see Graham, Knuth, and Patashnik, Chapter 7).

The series converges absolutely for

> abs(v) < 1/(exp(1)*(sqrt( (alpha^2+beta^2)/2 )));

>

Basic Facts about W

Applications

There are an extraordinary number of applications of W. This is not surprising, in retrospect: w*exp(w)=z is about the simplest transcendental equation that there is, that doesn't already have a solution in terms of known functions. A brief list of applications follows:

Stability analysis for delay differential equations; counting rooted, labelled trees; mathematical models of dozens of unrelated phenomena, including water movement in soil, Wien's displacement law (discussed in the new paper by S.R. Valluri, D.J. Jeffrey, and R.M. Corless, ``S ome Applications of the Lambert W function to Physics ''), combustion, epidemics, jet fuel consumption; and education.

Iteration of the exponential function (Euler, Carlsson, Hayes, Wright, ...)

The following problem is sometimes posed on math contests, and is answered on the FAQ for sci.math (or was, for a while).

> x^(x^(x^(`.`^(`.`^(`.`))))) = 2;

What is x? Once you've solved that one, now find x such that the left hand side is equal to 4. More generally, if we replace 2 by y, what is x as a function of y? What is y as a function of x? We are quickly led to the equation

> x^y = y;

This can be easily solved, in Maple:

> solve(%,y);

This tells us the limiting value of that iteration, when it converges. For more information, see the paper by Baker and Rippon in the Monthly from 1982.

Further, this iteration does not converge for all complex x. For some x, indeed, it settles to a two-cycle. For some, it settles to a three-cycle. Hayes first established that there were complex x-values giving rise to cycles in this iteration of arbitrary length. Drawing those x-values, in the right coordinate system, leads to a lovely image, where the regions of cycles are bordered by what looks like a fractal.

Complex Analysis

Branches

The Quadratrix of Hippias

We choose the branch cuts for W(k,z) by looking at the image of the negative real axis:

> restart;

> w := u + I*v;

> z := evalc(w*exp(w));

> x := evalc(Re(z));

> y := evalc(Im(z));

> solve(y=0,u);

> plot({seq([-v*cos(v)/sin(v),v,v=(k-1)*Pi*(1+(k-0.01)/1000)..k*Pi*(1-k/1000)],k=-3..4)},
view=[-4*Pi..4*Pi,-4*Pi..4*Pi],colour=black);

That is a subset of the curve known since antiquity as the Quadratrix of Hippias. The first curve ever named after its discoverer, Hippias of Elis (ca. 400BC), it can be used to square the circle, trisect the angle and even duplicate the cube. It is, therefore, not constructible by straightedge and compass.

But we only need those parts of it that satisfy x < 0. So we must look at the sign of x.

> X := simplify(subs(u=-v*cot(v),x));

> plot( X, v=-4*Pi..4*Pi, view=[-8..8,-2..2], numpoints=101, discont=true );

That allows us to deduce which pieces of the curve demark branch cuts for the Lambert W function. As a final filip,

> limit(X, v=0);

>

An Animation of a few values of W(k,z) for varying z and k (recondite, says David)

Now let's look at the complex branches of W. There is a symmetry relation amongst the branches, but it doesn't hold on the branch cuts unless we use a `signed zero'.

> restart;
alias(W=LambertW):

> fewvals := a -> [seq(evalf(W(k,a)),k=-6..5)];

> fewvals(-0.2):

> pts := a -> map( t->[evalc(Re(t)),evalc(Im(t))], fewvals(a));

> with(plots):

> N := 50;

> p := array(1..N);

> avals := [seq( evalf(-Pi/2*k/N),k=1..N)]:

> for i to N do p[i] := plot(pts(avals[i]),style=POINT,symbol=DIAMOND,symbolsize=30); od:

> plots[display]([seq(p[i],i=1..N)],insequence=true);

>

Riemann Surface

The following shows an animation of a section of the Riemann Surface for the Lambert W function. Strictly speaking, we have to prove that the technique works, by verifying that given (x,y,v) we can find u; that is, there is a bijection between a point on the 3d surface and the mapping z -> W(z) (i.e. (x,y) -> (u,v) ).

The Animated Riemann Surface for The Lambert W function

> restart;

> w := u+I*v;

> z := evalc(w*exp(w)):

> x := evalc(Re(z));

> y := evalc(Im(z));

> B := 6:

> f := theta -> plot3d([x,y,v], u=-6..1, v=-B..B, grid=[50,50], orientation=[theta,80], color=u):

> with(plots):

> bunch := [seq(f(10*k),k=-17..18)]:

> display(bunch, insequence=true, style=PATCHNOGRID, axes=NONE, view=[-1..1,-1..1,-B..B]);

>

>

A proof that the computed surface is correct

We must show that given (x,y,v) we can find u.

> restart;

> w := u+I*v;

> z := evalc(w*exp(w)):

> x := evalc(Re(z));

> y := evalc(Im(z));

> X/Y = normal( x/y );

> solve( %, u );

That solution is unique. Now we must investigate the condition y=0, which would have prevented us doing the division X/Y above.

> y=0;

One solution is v=0 (independently of u). If v=0, then u is determined by the x equation.

> X = eval( x, v=0);

There is a unique real solution if X >= 0, and two real solutions if -1/e <= X < 0, and no real solutions otherwise (this just follows from a graph).

Now, we should see if there are any possible solutions if v is not zero. If sin(v)=0, but v is not zero, then cos(v) = (-1)^k.

> 0 = eval(y, {sin(v)=0, cos(v)=(-1)^k});

This can happen only if u is -infinity.

Therefore our picture is good except on the line -1/e <= z < 0, when we have a spurious intersection.

>

Calculus

Derivatives

Differentiation of W(x)

> restart;
alias(W=LambertW);

> diff(W(x),x);

> normal(diff(W(x),x$5));

> diff(W(exp(x)),x);

> normal(diff(W(exp(x)),x$5));

> plot( {W(exp(x)),W(x)}, x =-1..2 );

>

Series

Now for some fun with series.

> restart;
alias(W=LambertW);

> series(W(x), x, 8);

The general term there is (-n)^(n-1)/n!. This is the exponential generating function for unrooted, labelled trees (apart from sign).

> add( (-n)^(n-1)/n!*x^n, n=1..7 );

Now a branch point series:

> series(W(x),x=-1/exp(1));

Let's get rid of those ugly square roots and exp(-1)'s:

> series(W((p^2/2-1)*exp(-1)), p);

Maple doesn't know the asymptotic series for W, at least not directly (because, to `asympt', W is an answer, not a question...)

> asympt(W(x),x);

But we can trick Maple into giving us what we want:

>

> alias( omega = RootOf( y + ln(y) - z , y ) );

> map( factor, asympt( omega, z ) );

That series is actually convergent (proof by de Bruijn). Comtet showed that the numbers that show up are Stirling numbers.

We will see later that omega is an interesting function itself; for now,

> solve( y + ln(y) = z, y );

Therefore, omega(z) = W(exp(z)), and so the asymptotics for W can be determined in Maple from the asymptotics of omega.

>

Integration

``Everybody knows'' that W(z) is not elementary. If it's true, where's the proof? Well, it hasn't been done, yet! James Davenport is currently working with David and me on completing the proof that W is not expressible in terms of elementary functions.

Nonetheless, we may integrate some expressions containing the Lambert W function, by the simple change of variable w = W(x) (and hence dx = (w+1)exp(w)dw ).

> restart;
alias(W=LambertW);

> int( sin(W(x)),x);

> plot(%,x=0..1);

> F := int( W(x)*ln(W(x)), x);

> eval( F, x=2.3 );

> plot( {evalc(Re(F)),evalc(Im(F))}, x=-1/10..1 , colour=black);

>

Computation

The method of choice in Maple to compute numerical values of W is Halley's method, which is a third-order variant of Newton's method. In fact, arbitrary-order methods for the computation of W exist, and are continually rediscovered;

> restart:

> f := w -> w*exp(w) - z;

> Halley := w -> w - f(w)/(D(f)(w) - 1/2*D(D(f))(w)*f(w)/D(f)(w));

> z := rand()/Float(1,12);

>

> Digits := 10;

> computed_w[0] := z;

> residual := computed_w[0]*exp(computed_w[0]) - z;

> Digits := 1;

>

> for k to 4 do Digits := 3*Digits; computed_w[k] := Halley( computed_w[k-1] ); od;

> Digits := Digits + 5;

> residual := computed_w[4]*exp(computed_w[4]) - z;

>

A surprising integral

Let us define the following strange-looking integral.

> restart:

Define the integrand, which depends on a parameter z.

> f := ((1-v*cot(v))^2 + v^2 )/(z + v*csc(v)*exp(-v*cot(v)));

Let's take a typical value of z, and plot the integrand:

> plot( eval(f,z=1), v=-Pi..Pi, -0.1..1.5 );

> F := Int( f, v=0..Pi )/Pi;

From R. William Gosper:

Finally, I must publicly "lambaste" Rob Corless for his preposterous allegation that

> W(z)/z = F;

Simply because Taylor expanding and numerically integrating gives

>

> series( F, z );

> evalf( eval( %, int=((m,n)->`evalf/int`(m,n,Digits,_Dexp))) );

> evalf( series( LambertW(z)/z, z, 7 ) );

>

...is no reason to go about endangering public morals and threatening massive unemployment via wild claims that implicit trascendental equations can be inverted with ordinary definite integrals. Sheesh!

Think of the numericists. Think of their families! (The numerical integrator converged slambangularly on these. This integral may well provide new approximations.)

RMC again: in fact, this integral is not the first integral expression for W. There is another, by Siewert and Burniston, which dates to 1975. Theirs is more complicated than the above, however.

After Bill Gosper's message, Steven Finch (at Mathsoft) pointed out that there is an integral in Bierens de Haan (1867) for this function (well, actually for T(x) = -W(-x))! It is expressed there in terms of infinite series, and does not use the notation W(x) or T(x). We are currently chasing this reference down.

Theorem : W(z)/z is a Stieltjes function , and hence its continued fraction has some well-understood properties (including the fact that all the poles of the subdiagonal Pade approximants of increasing degrees are real and interlace, for example). Proof: establish the integral above by use of the Cauchy integral formula and a transformation.

> alias(W=LambertW):

> series( W(z)/z, z, 8 );

> convert( %, ratpoly );

> fsolve( denom(%), z, complex );

> series( W(z)/z, z, 10 ):

> convert(%, ratpoly):

> fsolve( denom(%), z, complex );

Time permitting, we will explore the following formula for a bi-infinite family of zeros of the denominator of the integrand f:

> denom( f );

> zeros := (W(k,z)-W(m,z))/(2*I);

> subs( v=zeros, convert( denom(f), exp ) );

> simplify( % );

Petr Lisonek helped with the verification of that particular discovery.

>

The Wright function

The Wright function is defined by

> omega(z) = W( K(z), exp(z) );

where the unwinding number is defined by

> K := z -> ceil( (Im(z)-Pi)/(2*Pi) );

> plot(K(I*Pi*y),y=-4..4);

The reason that this is interesting is because

> ln(exp(z)) = z + 2*Pi*I*'K'(z);

Now we will see the relation of omega to the solution of the following equation:

> y + ln(y) = z;

It is clearly a cognate of W:

> _EnvAllSolutions := true;

> solve( y + ln(y) = z, y );

Maple says that all solutions of y + ln(y) = z are contained in that formula, for integers _NN1.

However, it turns out that not all choices for _NN1 lead to solutions.

Theorem : (Wright 1959, Siewert & Burniston 1975, Jeffrey, Corless & Hare 1995):
The solution of the equation y + ln y = z is

> y = piecewise( 'z=t+I*Pi and t <= -1', omega(z) and omega(z-2*Pi*I), 'z = t-I*Pi and t <= -1', "no solution", omega(z) );

That is, the solution of y + ln(y) = z is nearly unique, and can be expressed easily in terms of the single-valued function omega(z).

>

Applications

There are several applications of this function; indeed, some of them have been listed already as applications of W. But, because of superior numerical properties and other advantages, they are probably best listed as applications of directly. The first to come to our attention was the need in convex optimization. Consider the convex conjugate , namely, , of the function . Calculation shows that the conjugate is just , or . For large s, computing this function meant first exponentiation (making a REALLY large number, with possible overflow difficulties) and then taking (essentially) its logarithm. For some purposes, it is better to work directly with the Wright function.

%1. J.M. Borwein and A.S. Lewis, , CMS Advanced Books in Mathematics,.
@book{BorweinLewis:2000,
author = {Jonathan M.~Borwein and Adrian S.~Lewis},
title = {Convex Analysis and Nonlinear Optimization},
year = {2000},
publisher = {Canadian Mathematical Society},
series = {Advanced Books in Mathematics},
}

Another application is discussed in S.R. Valluri, D.J. Jeffrey, and R.M. Corless, ``Some Applications of the Lambert W function to Physics'', namely the conformal mapping for the fringing fields of a capacitor. We need to solve

> restart;

> 1 + z + exp(z) = zeta;

> solve(%, z );

This is better expressed (correctly accounting for branches, which Maple ignores above) as

> z = sort( zeta - 1 - omega( zeta - 1 ) );

>

Complex Analysis

Branches

Relations of the branches of W with .

We have that

> W(k,z) = omega( ln(z) + 2*Pi*I*k );

and by definition

> omega(z) = W(K(z),exp(z));

>

Riemann Surface

We examine the Riemann Surface for in the same manner that we looked at that for W.

Actual plot of the Riemann Surface for .

> restart:

> z := x+I*y; omega := mu + I*nu;

> x := evalc(Re(omega+ln(omega)));

> y := evalc(Im(omega+ln(omega)));

> plot3d( [x,y,mu], mu=-4..2, nu=-4..4, colour=black, axes=BOXED, style=PATCHNOGRID, labels=["x","y","mu"], view=[-2..1, -5..5, -5..3], grid=[200,200], style=POINT );

>

Proof that this computation works (incomplete)

We must check that given 3 of the coordinates, x, y, and mu as above, we can uniquely identify the nu, the fourth coordinate. If we can, then there is a bijection between the points on the 3-d surface and the mapping z to .

> restart:

> z := x+I*y; omega := mu + I*nu;

We ignore the case less than 0, which is the image of the branch cuts on z = t +/- I*Pi for t <= -1.

> x := evalc(Re(omega+ln(omega)));

> y := evalc(Im(omega+ln(omega)));

> solve( x=X, nu );

There are apparently two values of nu. However, if y and nu are not zero, (remember nu = 0 and mu < 0 corresponds to the image of the branch cut) then the equation with the arctan above uniquely determines the sign of nu, because y cannot be nu + arctan(nu,mu) and -nu + arctan(-nu,mu) simultaneously. If nu is zero, then

> Y = 0 + arctan(0,mu);

This means that y must be zero, if mu > 0, or Pi, if mu < 0. The case when y=+/-Pi must be handled separately (it is in the branch cut case that we are ignoring for the moment). If nu = 0 and mu > 0, the x-equation becomes

> _EnvAllSolutions := true;

> X = mu + ln(mu);

Since mu and x are real, we get something that we can just plot:

> plot([mu+ln(mu),mu,mu=0.0001..3]);

Now it only remains to examine the branch cut, and the case y=0 but nu not zero.

> penultimate_case := eval([X=x,0=y]);

If there are any zeros of this at all, then there are two values of nu corresponding to every value of mu, because these equations are symmetric under the change of sign nu -> -nu; this will produce a double root at nu=0, of course, but this is excluded.

A short inspection of the sign of the tangent in each quadrant convinces us that for real nu and mu, there is no solution to

> 0 = arctan(nu,mu) + nu;

unless nu = 0.

At long last, considering the branches, we have

> restart;

> omega(Z) = piecewise(Z=t+I*Pi, W[0](-exp(t)),Z=t-I*Pi,W[-1](-exp(t)));

and in both cases we know that nu = 0. Therefore, knowing x, y, and mu in this case means we also know nu; therefore each point on these lines corresponds uniquely to a piece of the mapping from z to omega, and therefore is a valid representation of the Riemann Surface for omega.

>

>

Calculus

Derivatives

The derivative of is slightly simpler than that of W.

> restart;
alias(omega=RootOf(y+ln(y)=z,y));

> normal( diff( omega, z ) );

> normal( diff( omega, z$5 ) );

Those coefficients are 2nd order Eulerian numbers. Don Knuth was the one to recognize this fact.

>

Series

Series for are likewise simpler.

> series( omega, z=1 );

> alias( w=LambertW(exp(a)) );

> series( omega, z=a );

> series( omega, z=-1+I*Pi );

Error, (in series/RootOf) unable to compute series

So, Maple doesn't know how to do that series in this form. In fact, that series is completely known, and is related to Stirling's formula for n! (see Corless, Jeffrey, and Knuth 1997).

> alias(W=LambertW);

> series(W(-exp(-1-p^2/2)),p,10);

> asympt(GAMMA(n),n);

In his recent talk at Waterloo, Don Knuth shows a connection between the even coefficients of the W series with Ramanujan's function.

>

Integration

Again, we can evaluate some integrals containing the Wright function, by a change of variables. Maple has just been taught this change of variables.

> restart;
read "wright.mpl";
with(Wright);

> F := Int( sin(omega(z)), z );

> value(F);

Laplace Transform of the Wright function

> with(inttrans):

> laplace( omega(z), z, s );

> F := Int( omega(z)*exp(-s*z), z );

> Laplace_integral_omega := simplify(value(expand(simplify(F))));

I don't know how to simplify that any further with Maple! However, the answer is really pretty simple---I leave it to you to work out for yourself.

>

Computation

As before, we have a family of arbitrary-order numerical iterative methods to compute the Wright function.

> restart;
alias( omega=RootOf(y+ln(y)-z,y),
w = LambertW(exp(a)) );

> series( omega, z=a, 4 );

> fourth := (w,z) -> (w + w/(1+w)*(z-w-ln(w)) + w/(1+w)^3*(z-w-ln(w))^2/2 - w*(2*w-1)/(1+w)^5*(z-w-ln(w))^3/6);

> z := rand()/Float(1,12);

> computed_omega[0] := z;

> Digits := 1;

> for i to 3 do
Digits := 4*Digits;
computed_omega[i] := fourth( computed_omega[i-1], z );
od;

> Digits := Digits + 5;

> z - computed_omega[3] - ln(computed_omega[3]);

Near the discontinuities (remember the Riemann surface) more care has to be taken; we have found a regularizing transformation that allows us to accurately iterate in this fashion.

> restart;

> read "wright.mpl";

> with(Wright):

> z := .5584587190;

> Digits := 400;

> omega(z);

> omega(z) + ln(omega(z)) - z;

>

Concluding Remarks

Thank you for listening to me obsess ; I hope that you have enjoyed some of these results, and I also hope that you will feel, when you see W (or omega) in a solution to one of your problems, that you have an actual answer.