{"id":754,"date":"2013-09-02T12:00:33","date_gmt":"2013-09-02T17:00:33","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=754"},"modified":"2017-02-24T20:16:32","modified_gmt":"2017-02-25T01:16:32","slug":"the-lambert-w-function","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/02\/the-lambert-w-function\/","title":{"rendered":"The Lambert W Function"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p>The Lambert W function deserves to be better known.  It pops up in all sorts of places.   And our MATLAB function for evaluating the function is a beautiful use of the Halley method.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#bfba4e2d-e049-45a6-8285-fe8b51d69ce7\">Johann Lambert<\/a><\/li><li><a href=\"#92be7f49-9302-41e4-876a-d0c41789f465\">The Maple Connection<\/a><\/li><li><a href=\"#a2d91197-cb60-4188-b7f0-36dad17ae29d\">An elementary function<\/a><\/li><li><a href=\"#5ae6f42f-d5da-456f-8011-201783a5232c\">Functional inverse<\/a><\/li><li><a href=\"#0294079d-9227-4a93-846e-1773193df571\">Halley's method<\/a><\/li><li><a href=\"#3c285ad4-958b-4841-8dc8-15df5b9a37a7\">Application to Lambert W<\/a><\/li><li><a href=\"#d95f31b1-688c-4a01-9759-33ae9481134e\">Starting values<\/a><\/li><li><a href=\"#4d21d68a-439f-4838-893e-8864a00c3adf\">Stopping values<\/a><\/li><li><a href=\"#c1ea75cb-227b-4aec-959b-6d5b66c07d54\">Lambert_W.m<\/a><\/li><li><a href=\"#3ffe1fd2-6ec1-44aa-8eaa-a389397d8bd0\">Plot Lambert W<\/a><\/li><li><a href=\"#94a94b80-116f-4c41-9716-f7a6c8268310\">Acknowledgements<\/a><\/li><li><a href=\"#01da2f83-4952-4750-8cc1-632969d9e5f6\">Lambert W Poster<\/a><\/li><li><a href=\"#ef91b7df-a138-4b5d-8abf-c8ee368c3ef1\">Further Reading<\/a><\/li><\/ul><\/div><h4>Johann Lambert<a name=\"bfba4e2d-e049-45a6-8285-fe8b51d69ce7\"><\/a><\/h4><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/JHLambert.jpg\" alt=\"\"> <\/p><p>Johann Henrich Lambert was born in 1728 in Mulhouse, which was then in Switzerland, but is now in France. He was a contemporary of the famous Swiss mathematician Leonard Euler.  He had incredibly wide ranging interests. He was the first to prove that $\\pi$ was irrational.  He introduced the hyperbolic trig functions.  He made such important contributions to optics that a photometric unit for luminance is named in his honor.  He also made contributions in physics, astronomy, logic, and philosophy.<\/p><p>Lambert did not explicitly describe the function that is now named after him.  He considered various equations, including ones of the form $x = q + x^m$.  Euler followed up with papers about Lambert's work. The ideas involved in analyzing such equations underlie what we today call the Lambert W function.<\/p><h4>The Maple Connection<a name=\"92be7f49-9302-41e4-876a-d0c41789f465\"><\/a><\/h4><p>In the early 1990s developers of the Maple symbolic algebra software, including Rob Corless, Gaston Gonnet, D.E.G. Hare, and David Jeffrey, came to realize that this function had been rediscovered many times over the years.  They christened the function \"W\", published a note in the Maple Technical Newsletter, and began to develop a bibliography. Don Knuth joined the group by contributing applications involving enumerations of trees. In 1996 the five of them wrote a 32 page paper that has become the definitive reference on the Lambert W function.  Their bibliography contains 72 entries.<\/p><h4>An elementary function<a name=\"a2d91197-cb60-4188-b7f0-36dad17ae29d\"><\/a><\/h4><p>Let's plot an elementary function, using $w$ as the independent variable. The function is $w \\mbox{e}^w$.<\/p><pre class=\"codeinput\">   f = @(w) w.*exp(w);\r\n   ezplot(f,[-4 1])\r\n   xlabel(<span class=\"string\">'w'<\/span>)\r\n   axis([-4 1 -.5 1])\r\n   axis <span class=\"string\">square<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lambertw_blog_01.png\" alt=\"\"> <p>You can see the function reaches a global minimum at $w = -1$.  The value of this minimum plays an important role in this story, so let's give it a name.<\/p><p>$$ \\bar{\\mbox{e}} = -1\/\\mbox{e} $$<\/p><pre class=\"codeinput\">   ebar = -1\/exp(1)\r\n   line([-1,-1],[ebar,ebar],<span class=\"string\">'marker'<\/span>,<span class=\"string\">'.'<\/span>,<span class=\"string\">'markersize'<\/span>,18,<span class=\"string\">'color'<\/span>,<span class=\"string\">'black'<\/span>)\r\n   set(gca,<span class=\"string\">'ytick'<\/span>,[ebar,0,0.5,1])\r\n   set(gca,<span class=\"string\">'yticklabel'<\/span>,{<span class=\"string\">'-1\/e'<\/span>,<span class=\"string\">'0'<\/span>,<span class=\"string\">'0.5'<\/span>,<span class=\"string\">'1'<\/span>})\r\n<\/pre><pre class=\"codeoutput\">ebar =\r\n  -0.367879441171442\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lambertw_blog_02.png\" alt=\"\"> <h4>Functional inverse<a name=\"5ae6f42f-d5da-456f-8011-201783a5232c\"><\/a><\/h4><p>The Lambert W function, $W(x)$, is the functional inverse of $w \\mbox{e}^w$. In other words, $w$ = W($x$) is the solution of the equation<\/p><p>$$ x = w \\mbox{e}^w $$<\/p><p>The graph of $W(x)$ is conceptually obtained by labeling the vertical axis of the above graph with an <tt>'x'<\/tt> and then interchanging the horizontal and vertical axes.  (You can sneak a peak at the actual plot below, generated after we see how to compute the function.)<\/p><p>The point $x = \\bar{\\mbox{e}} = -1\/\\mbox{e}$ is crucial. On the interval $\\bar{\\mbox{e}} &lt;= x &lt;= 0$ there are two real solutions and so $W(x)$ is double valued. One solution is increasing, is defined for all $x &gt;= \\bar{\\mbox{e}}$, is known is the <i>primary branch<\/i>, and is denoted by $W_0(x)$. The other solution is decreasing, is defined only for negative $x$, has a pole at 0, and is denoted by $W_{-1}(x)$. (There are other, complex-valued, branches, $W_k(x)$, which we are not considering here.)<\/p><h4>Halley's method<a name=\"0294079d-9227-4a93-846e-1773193df571\"><\/a><\/h4><p>Halley's method is a method for finding a zero of a real-valued function with a continuous and easily computed second derivative. It is named after the British astronomer who is better known for discovering a comet.<\/p><p>Newton's method approximates a function locally by a linear function with the same slope and steps to the zero of that approximation. Halley's method approximates the function locally by a hyperbola with the same slope and curvature and steps to the nearest zero of that approximation. The iteration is<\/p><p>$$w_{n+1}=w_n-f(w_n)\/\\left[f'(w_n)-\\frac{f(w_n)f''(w_n)}{2f'(w_n)}\\right]$$<\/p><p>You can see Newton's method by ignoring the second term in the square brackets.  You can see that Halley's method is particularly efficient if the ratio $f''(w)\/f'(w)$ happens to take a simple form.  But we do not regard Halley's method as a general purpose method because it does require knowledge of the second derivative.<\/p><h4>Application to Lambert W<a name=\"3c285ad4-958b-4841-8dc8-15df5b9a37a7\"><\/a><\/h4><p>Computing $W(x)$ for a given value of $x$ requires solving the equation<\/p><p>$$ f(w) = w \\mbox{e}^w - x = 0 $$<\/p><p>for $w$.  We can use Halley's method because the required derivatives are readily available.<\/p><p>$$ f'(w) = w \\mbox{e}^w + \\mbox{e}^w = \\mbox{e}^w (w + 1) $$<\/p><p>$$ f''(w) = w \\mbox{e}^w + w \\mbox{e}^w + \\mbox{e}^w = \\mbox{e}^w (w + 2) $$<\/p><p>Note that<\/p><p>$$ f''(w)\/f'(w) = (w+2)\/(w+1) $$<\/p><p>This gives us the following iterative step for computing Lambert W at a given value of $x$.<\/p><p>$$ f_n = w_n \\mbox{e}^{w_n} - x $$<\/p><p>$$w_{n+1}=w_n-f_n\/\\left[\\mbox{e}^{w_n}(w_n+1)-(w_n+2)f_n\/(2 w_n+2)\\right]$$<\/p><h4>Starting values<a name=\"d95f31b1-688c-4a01-9759-33ae9481134e\"><\/a><\/h4><p>Various people have written programs before me to compute Lambert W using Halley's method.  They have paid quite a bit of attention to starting values that reduce the eventual number of iterations. I have decided to take a different approach.  I am willing to have the computation take a few more iterations if the program can be simplified. In fact, I use the simplest possible initial values.<\/p><p>For the primary branch, start with $w = 1$ for any $x$.  Do not try to approximate $W(x)$ at all.  The iteration turns out to converge to double precision floating point accuracy in about a half dozen iterations for any $x$, except values very near the branch point at $\\bar{\\mbox{e}}$. More accurate starting values could cut the number of iterations in half, but I prefer the simpler code.<\/p><p>For the lower branch, start with any value less than $-1$.  I have chosen $w = -2$, but I don't think it matters very much. Again, about a half dozen iterations are required unless $x$ is very close to the branch point at $\\bar{\\mbox{e}}$ or the pole at $0$.<\/p><h4>Stopping values<a name=\"4d21d68a-439f-4838-893e-8864a00c3adf\"><\/a><\/h4><p>Halley's method is cubically convergent near a simple root, which this is. Once we're close, we get closer very fast.  If two successive iterates agree to <tt>1.0e-8<\/tt>, then you can be pretty sure the second one is accurate to <tt>1.0e-16<\/tt>.  I could probably even replace <tt>1.e-8<\/tt> by <tt>eps^(1\/3) = 6e-6<\/tt>.<\/p><h4>Lambert_W.m<a name=\"c1ea75cb-227b-4aec-959b-6d5b66c07d54\"><\/a><\/h4><p>Here is the entire M-file for <tt>function Lambert_W<\/tt>.  This will be available from the MATLAB Central File Exchange in due course. (Added 9\/16\/2013: See <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/43419-the-lambert-w-function\">this link<\/a>)<\/p><pre class=\"codeinput\">   type <span class=\"string\">Lambert_W<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction w = Lambert_W(x,branch)\r\n% Lambert_W  Functional inverse of x = w*exp(w).\r\n% w = Lambert_W(x), same as Lambert_W(x,0)\r\n% w = Lambert_W(x,0)  Primary or upper branch, W_0(x)\r\n% w = Lambert_W(x,-1)  Lower branch, W_{-1}(x)\r\n%\r\n% See: https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/02\/the-lambert-w-function\/\r\n\r\n% Copyright 2013 The MathWorks, Inc.\r\n\r\n% Effective starting guess\r\nif nargin &lt; 2 || branch ~= -1\r\n   w = ones(size(x));  % Start above -1\r\nelse  \r\n   w = -2*ones(size(x));  % Start below -1\r\nend\r\nv = inf*w;\r\n\r\n% Haley's method\r\nwhile any(abs(w - v).\/abs(w) &gt; 1.e-8)\r\n   v = w;\r\n   e = exp(w);\r\n   f = w.*e - x;  % Iterate to make this quantity zero\r\n   w = w - f.\/((e.*(w+1) - (w+2).*f.\/(2*w+2)));\r\nend\r\n\r\n<\/pre><h4>Plot Lambert W<a name=\"3ffe1fd2-6ec1-44aa-8eaa-a389397d8bd0\"><\/a><\/h4><p>We can now use this code to produce a plot of both branches of the Lambert W function.<\/p><pre class=\"codeinput\"><span class=\"comment\">% Primary branch<\/span>\r\n\r\n   x = ebar:.01:1;\r\n   plot(x,Lambert_W(x,0))\r\n\r\n<span class=\"comment\">% Lower branch<\/span>\r\n\r\n   hold <span class=\"string\">on<\/span>\r\n   x = ebar:.01:0;\r\n   plot(x,Lambert_W(x,-1),<span class=\"string\">'r'<\/span>)\r\n   hold <span class=\"string\">off<\/span>\r\n\r\n<span class=\"comment\">% Annotation<\/span>\r\n\r\n   axis([-.5 1 -4 1])\r\n   axis <span class=\"string\">square<\/span>\r\n   line([ebar ebar],[-1 -1],<span class=\"string\">'marker'<\/span>,<span class=\"string\">'.'<\/span>,<span class=\"string\">'markersize'<\/span>,18,<span class=\"string\">'color'<\/span>,<span class=\"string\">'black'<\/span>)\r\n   set(gca,<span class=\"string\">'xtick'<\/span>,[ebar 0 .5 1])\r\n   set(gca,<span class=\"string\">'xticklabel'<\/span>,{<span class=\"string\">'-1\/e'<\/span>,<span class=\"string\">'0'<\/span>,<span class=\"string\">'0.5'<\/span>,<span class=\"string\">'1'<\/span>})\r\n   set(gca,<span class=\"string\">'ytick'<\/span>,-4:1)\r\n   xlabel(<span class=\"string\">'x'<\/span>)\r\n   ylabel(<span class=\"string\">'w'<\/span>)\r\n   title(<span class=\"string\">'Lambert W'<\/span>)\r\n   legend({<span class=\"string\">'W_0'<\/span>,<span class=\"string\">'W_{-1}'<\/span>},<span class=\"string\">'location'<\/span>,<span class=\"string\">'southeast'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lambertw_blog_03.png\" alt=\"\"> <h4>Acknowledgements<a name=\"94a94b80-116f-4c41-9716-f7a6c8268310\"><\/a><\/h4><p>Thanks to Peter Costa for reawakening my interest and to Rob Corless for keeping me up to date as always.<\/p><h4>Lambert W Poster<a name=\"01da2f83-4952-4750-8cc1-632969d9e5f6\"><\/a><\/h4><p>Follow this link to a poster about the Lambert W function produced by the group at the University of Western Ontario led by Corless and Jeffrey.<\/p><h4>Further Reading<a name=\"ef91b7df-a138-4b5d-8abf-c8ee368c3ef1\"><\/a><\/h4><p>Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, Donald (1996). \"On the Lambert W function\". Advances in Computational Mathematics (Berlin, New York: Springer-Verlag) 5: 329-359. <a href=\"http:\/\/link.springer.com\/article\/10.1007%2FBF02124750\">&lt;http:\/\/link.springer.com\/article\/10.1007%2FBF02124750<\/a>&gt;<\/p><p>Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, Donald (1996). \"On the Lambert W function\". PDF Preprint. <a href=\"https:\/\/cs.uwaterloo.ca\/research\/tr\/1993\/03\/W.pdf\">https:\/\/cs.uwaterloo.ca\/research\/tr\/1993\/03\/W.pdf<\/a><\/p><p>Wikipedia article. <a href=\"http:\/\/en.wikipedia.org\/wiki\/Lambert_W_function\">&lt;http:\/\/en.wikipedia.org\/wiki\/Lambert_W_function<\/a>&gt;<\/p><p>NIST Digital Library of Mathematical Functions. <a href=\"http:\/\/dlmf.nist.gov\/4.13\">&lt;http:\/\/dlmf.nist.gov\/4.13<\/a>&gt;<\/p><p>Halley's Method: Dahlquist, G., Bjorck, A. <i>Numerical Methods in Scientific Computing<\/i>, Volume I, page 648, SIAM, Philadelphia, 2008.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_e4c4463a660d4d20a71db185a329c2fe() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='e4c4463a660d4d20a71db185a329c2fe ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' e4c4463a660d4d20a71db185a329c2fe';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2013 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_e4c4463a660d4d20a71db185a329c2fe()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2013b<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; R2013b<br><\/p><\/div><!--\r\ne4c4463a660d4d20a71db185a329c2fe ##### SOURCE BEGIN #####\r\n%% The Lambert W Function\r\n% The Lambert W function deserves to be better known.  It pops up in all\r\n% sorts of places.   And our MATLAB function for evaluating the function\r\n% is a beautiful use of the Halley method.\r\n\r\n%% Johann Lambert\r\n%\r\n% <<JHLambert.jpg>>\r\n%\r\n% Johann Henrich Lambert was born in 1728 in Mulhouse, which was then in\r\n% Switzerland, but is now in France. He was a contemporary of the famous\r\n% Swiss mathematician Leonard Euler.  He had incredibly wide ranging interests.\r\n% He was the first to prove that $\\pi$ was irrational.  He introduced the\r\n% hyperbolic trig functions.  He made such important contributions to optics\r\n% that a photometric unit for luminance is named in his honor.  He also made\r\n% contributions in physics, astronomy, logic, and philosophy.\r\n\r\n%%\r\n% Lambert did not explicitly describe the function that is now named after\r\n% him.  He considered various equations, including ones of the form\r\n% $x = q + x^m$.  Euler followed up with papers about Lambert's work.\r\n% The ideas involved in analyzing such equations underlie what we today\r\n% call the Lambert W function.\r\n\r\n%% The Maple Connection\r\n% In the early 1990s developers of the Maple symbolic algebra software,\r\n% including Rob Corless, Gaston Gonnet, D.E.G. Hare, and David Jeffrey,  \r\n% came to realize that this function had been rediscovered many times\r\n% over the years.  They christened the function \"W\", published a note\r\n% in the Maple Technical Newsletter, and began to develop a bibliography.\r\n% Don Knuth joined the group by contributing applications involving\r\n% enumerations of trees.\r\n% In 1996 the five of them wrote a 32 page paper that has become the definitive\r\n% reference on the Lambert W function.  Their bibliography contains 72 entries.\r\n\r\n%% An elementary function\r\n% Let's plot an elementary function, using $w$ as the independent variable.\r\n% The function is $w \\mbox{e}^w$.\r\n%\r\n   f = @(w) w.*exp(w);\r\n   ezplot(f,[-4 1])\r\n   xlabel('w')\r\n   axis([-4 1 -.5 1])\r\n   axis square\r\n\r\n%%\r\n% You can see the function reaches a global minimum at $w = -1$.  The value\r\n% of this minimum plays an important role in this story, so let's give it\r\n% a name.\r\n%\r\n% $$ \\bar{\\mbox{e}} = -1\/\\mbox{e} $$\r\n%\r\n   ebar = -1\/exp(1)\r\n   line([-1,-1],[ebar,ebar],'marker','.','markersize',18,'color','black')\r\n   set(gca,'ytick',[ebar,0,0.5,1])\r\n   set(gca,'yticklabel',{'-1\/e','0','0.5','1'})\r\n\r\n%% Functional inverse\r\n% The Lambert W function, $W(x)$, is the functional inverse of $w \\mbox{e}^w$.\r\n% In other words, $w$ = W($x$) is the solution of the equation\r\n%\r\n% $$ x = w \\mbox{e}^w $$\r\n%\r\n% The graph of $W(x)$ is conceptually obtained by labeling the\r\n% vertical axis of the above graph with an |'x'| and then interchanging the\r\n% horizontal and vertical axes.  (You can sneak a peak at the actual plot\r\n% below, generated after we see how to compute the function.)\r\n\r\n%%\r\n% The point $x = \\bar{\\mbox{e}} = -1\/\\mbox{e}$ is crucial.\r\n% On the interval $\\bar{\\mbox{e}} <= x <= 0$ there are two real solutions\r\n% and so $W(x)$ is double valued.\r\n% One solution is increasing, is defined for all $x >= \\bar{\\mbox{e}}$,\r\n% is known is the _primary branch_, and is denoted by $W_0(x)$. \r\n% The other solution is decreasing, is defined only for negative $x$,\r\n% has a pole at 0, and is denoted by $W_{-1}(x)$. \r\n% (There are other, complex-valued, branches, $W_k(x)$, which we are not\r\n% considering here.)\r\n\r\n%% Halley's method\r\n% Halley's method is a method for finding a zero of a real-valued\r\n% function with a continuous and easily computed second derivative.\r\n% It is named after the British astronomer who is better known for\r\n% discovering a comet.\r\n\r\n%%\r\n% Newton's method approximates a function locally by a linear function\r\n% with the same slope and steps to the zero of that approximation.\r\n% Halley's method approximates the function locally by a hyperbola with the\r\n% same slope and curvature and steps to the nearest zero of that approximation.\r\n% The iteration is\r\n%\r\n% $$w_{n+1}=w_n-f(w_n)\/\\left[f'(w_n)-\\frac{f(w_n)f''(w_n)}{2f'(w_n)}\\right]$$\r\n%\r\n\r\n%%\r\n% You can see Newton's method by ignoring the second term in the square\r\n% brackets.  You can see that Halley's method is particularly efficient\r\n% if the ratio $f''(w)\/f'(w)$ happens to take a simple form.  But we do not\r\n% regard Halley's method as a general purpose method because it does\r\n% require knowledge of the second derivative.\r\n\r\n%% Application to Lambert W\r\n% Computing $W(x)$ for a given value of $x$ requires solving\r\n% the equation\r\n%\r\n% $$ f(w) = w \\mbox{e}^w - x = 0 $$\r\n%\r\n% for $w$.  We can use Halley's method because the required derivatives\r\n% are readily available.\r\n%\r\n% $$ f'(w) = w \\mbox{e}^w + \\mbox{e}^w = \\mbox{e}^w (w + 1) $$\r\n%\r\n%\r\n% $$ f''(w) = w \\mbox{e}^w + w \\mbox{e}^w + \\mbox{e}^w = \\mbox{e}^w (w + 2) $$\r\n%\r\n\r\n%%\r\n% Note that\r\n%\r\n% $$ f''(w)\/f'(w) = (w+2)\/(w+1) $$\r\n%\r\n\r\n%%\r\n% This gives us the following iterative step for computing Lambert W at\r\n% a given value of $x$.\r\n%\r\n% $$ f_n = w_n \\mbox{e}^{w_n} - x $$\r\n%\r\n% $$w_{n+1}=w_n-f_n\/\\left[\\mbox{e}^{w_n}(w_n+1)-(w_n+2)f_n\/(2 w_n+2)\\right]$$\r\n%\r\n\r\n%% Starting values\r\n% Various people have written programs before me to compute Lambert W using\r\n% Halley's method.  They have paid quite a bit of attention to starting\r\n% values that reduce the eventual number of iterations.\r\n% I have decided to take a different approach.  I am willing to have the\r\n% computation take a few more iterations if the program can be simplified.\r\n% In fact, I use the simplest possible initial values.\r\n\r\n%% \r\n% For the primary branch, start with $w = 1$ for any $x$.  Do not try to\r\n% approximate $W(x)$ at all.  The iteration turns out to converge to\r\n% double precision floating point accuracy in about a half dozen iterations\r\n% for any $x$, except values very near the branch point at $\\bar{\\mbox{e}}$.\r\n% More accurate starting values could cut the number of iterations in half,\r\n% but I prefer the simpler code.\r\n\r\n%%\r\n% For the lower branch, start with any value less than $-1$.  I have chosen\r\n% $w = -2$, but I don't think it matters very much.\r\n% Again, about a half dozen iterations are required unless $x$ is very close\r\n% to the branch point at $\\bar{\\mbox{e}}$ or the pole at $0$.\r\n\r\n%% Stopping values\r\n% Halley's method is cubically convergent near a simple root, which this is.\r\n% Once we're close, we get closer very fast.  If two successive iterates\r\n% agree to |1.0e-8|, then you can be pretty sure the second one is accurate\r\n% to |1.0e-16|.  I could probably even replace |1.e-8| by |eps^(1\/3) = 6e-6|.\r\n\r\n%% Lambert_W.m\r\n% Here is the entire M-file for |function Lambert_W|.  This will be \r\n% available from the MATLAB Central File Exchange in due course.\r\n% (Added 9\/16\/2013: See\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/43419-the-lambert-w-function\r\n% this link>)\r\n\r\n   type Lambert_W\r\n\r\n%% Plot Lambert W\r\n% We can now use this code to produce a plot of both branches of the\r\n% Lambert W function.\r\n\r\n% Primary branch\r\n\r\n   x = ebar:.01:1;\r\n   plot(x,Lambert_W(x,0))\r\n\r\n% Lower branch\r\n\r\n   hold on\r\n   x = ebar:.01:0;\r\n   plot(x,Lambert_W(x,-1),'r')\r\n   hold off\r\n\r\n% Annotation\r\n\r\n   axis([-.5 1 -4 1])\r\n   axis square\r\n   line([ebar ebar],[-1 -1],'marker','.','markersize',18,'color','black')\r\n   set(gca,'xtick',[ebar 0 .5 1])\r\n   set(gca,'xticklabel',{'-1\/e','0','0.5','1'})\r\n   set(gca,'ytick',-4:1)\r\n   xlabel('x')\r\n   ylabel('w')\r\n   title('Lambert W')\r\n   legend({'W_0','W_{-1}'},'location','southeast')\r\n\r\n%% Acknowledgements\r\n% Thanks to Peter Costa for reawakening my interest and to Rob Corless for\r\n% keeping me up to date as always.\r\n\r\n%% Lambert W Poster\r\n% Follow <http:\/\/www.orcca.on.ca\/LambertW\/ this link> to a poster about\r\n% the Lambert W function produced by the group at the University of\r\n% Western Ontario led by Corless and Jeffrey.\r\n\r\n%% Further Reading\r\n% Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, Donald (1996).\r\n% \"On the Lambert W function\". Advances in Computational Mathematics\r\n% (Berlin, New York: Springer-Verlag) 5: 329-359.\r\n% <http:\/\/link.springer.com\/article\/10.1007%2FBF02124750 \r\n% http:\/\/link.springer.com\/article\/10.1007%2FBF02124750>\r\n\r\n%%\r\n% Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, Donald (1996).\r\n% \"On the Lambert W function\". PDF Preprint.\r\n% <https:\/\/cs.uwaterloo.ca\/research\/tr\/1993\/03\/W.pdf\r\n% https:\/\/cs.uwaterloo.ca\/research\/tr\/1993\/03\/W.pdf>\r\n\r\n%%\r\n% Wikipedia article.\r\n% <http:\/\/en.wikipedia.org\/wiki\/Lambert_W_function\r\n% http:\/\/en.wikipedia.org\/wiki\/Lambert_W_function>\r\n\r\n%%\r\n% NIST Digital Library of Mathematical Functions.\r\n% <http:\/\/dlmf.nist.gov\/4.13 http:\/\/dlmf.nist.gov\/4.13>\r\n\r\n%%\r\n% Halley's Method: Dahlquist, G., Bjorck, A.\r\n% _Numerical Methods in Scientific Computing_, Volume I, page 648,\r\n% SIAM, Philadelphia, 2008. \r\n\r\n##### SOURCE END ##### e4c4463a660d4d20a71db185a329c2fe\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/JHLambert.jpg\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>The Lambert W function deserves to be better known.  It pops up in all sorts of places.   And our MATLAB function for evaluating the function is a beautiful use of the Halley method.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/02\/the-lambert-w-function\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,4,16,17],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/754"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=754"}],"version-history":[{"count":12,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/754\/revisions"}],"predecessor-version":[{"id":2360,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/754\/revisions\/2360"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=754"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=754"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=754"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}