I began working on Lambert's W a while back, as I needed it to numerically
evaluate something. So far, I have a crude numerical implementation, not yet
ready for prime time, but it might work in some cases:
LambertW(x):=if not numberp(x) then nounify(LambertW)(x)
else if x<=float(-1/%e) then minf
else if x<=-0.3 then float(-1.666666667+2.331643981*(x+.3678794412)\
^(1/2)-1.812187885*x+1.936631115*(x+.3678794412)^(3/2)-2.353551201*\
(x+.3678794412)^2+3.066858901*(x+.3678794412)^(5/2)-4.175335600*(x+\
.3678794412)^3+5.858023730*(x+.3678794412)^(7/2)-8.401032216*(x+\
.3678794412)^4+12.25075350*(x+.3678794412)^(9/2)-18.10069700*(x+\
.3678794412)^5+27.02904480*(x+.3678794412)^(11/2))
else if x<=float(1/%e) then float(sum((-i)^(i-1)/i!*x^i,i,1,3*fpprec))
else block
(
[w,we,w1e],
w:if x<=500.0 then block([lx1:log(x+1.0)], 0.665*(1+0.0195*lx1)*lx1+0.04)
else log(x-4.0) - (1.0 - 1.0/log(x))*log(log(x)),
for i thru 100 do
(
we:w * exp(w),
w1e:(w + 1) * exp(w),
if (x - we) / w1e < 10.0^(1-fpprec) then i:100
else w:w - (we - x) / (w1e - (w+2) * (we-x) / (2*w+2))
),
w
)$
The derivative of Lambert's W can be defined as
gradef(LambertW(x),LambertW(x)/(1+LambertW(x))/x);
I've been meaning to spend some time on this and turn it into something
that's more "production quality".
Viktor
-----Original Message-----
From: maxima-bounces at math.utexas.edu [mailto:maxima-bounces at math.utexas.edu]
On Behalf Of Oliver Kullmann
Sent: Tuesday, May 20, 2008 12:41 PM
To: maxima at math.utexas.edu
Subject: Lambert W function?
Hi,
it seems Maxima doesn't know the Lambert W function, see
http://en.wikipedia.org/wiki/Lambert_W_function
With that function one could enhance "solve" to handle
also for example e^x - x = a.
(I needed such a case, and my old Mupad implementation
returned the solution using that W-function, which I found
helpful.)
Oliver
_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima