Subject: Integrator loops endlessly for lambert_w(1/x)
From: Viktor T. Toth
Date: Fri, 26 Dec 2008 09:32:57 -0500
A while back, I proposed using a different expression for the derivative of
Lambert's W, as it leads to better simplification:
gradef(lambert_w(x),lambert_w(x)/(1+lambert_w(x))/x);
The downside of using this version of the derivative is that it's undefined
at x=0, and (for now, anyway) Maxima cannot compute its limit either. On the
other hand, it leads to expressions that can be simplified much more
readily, due to the fact that lambert_w doesn't appear in exponents.
If I use this definition as the derivative of Lambert's W,
integrate(lambert_w(x),x) just gives back the noun form, as Maxima is unable
to compute the integral, but no stack overflow occurs.
The integral itself should be
x*lambert_w(1/x)+expintegral_e(1,lambert_w(1/x));
which can be verified, for instance, by
diff(%,x);
ratsubst(lambert_w(1/x)*x,exp(-lambert_w(1/x)),%);
Hmmm, perhaps there is a way to teach the simplifier that z =
lambert_w(z)*exp(lambert_w(z))?
Viktor
-----Original Message-----
From: maxima-bounces at math.utexas.edu [mailto:maxima-bounces at math.utexas.edu]
On Behalf Of Dieter Kaiser
Sent: Thursday, December 25, 2008 11:30 PM
To: maxima at math.utexas.edu
Subject: Integrator loops endlessly for lambert_w(1/x)
Today I have recognized that David Billinghurst has added a test to
rtest_lambert_w.mac which fails:
integrate(lambert_w(1/x),x)
The integrator loops endlessly.
I had a long search for the problem. The integrator tries to do an
integration-by-parts and generates the following integrand:
%e^(-lambert_w(1/x))/(x*(1+lambert_w(1/x)))
This is the expression we get when we type it in at the Maxima prompt
(The numbers enumerate the different parts of the expression):
((MTIMES SIMP)
1* ((MEXPT SIMP) $X -1)
2* ((MEXPT SIMP) $%E
((MTIMES SIMP) -1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))))
3* ((MEXPT SIMP)
((MPLUS SIMP) 1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))) -1))
Then I have looked into the integrator. In the integrator the expression
is a bit different ordered. The expression is generated in the routine
partial-integration and transfered to the integrator. This is the
expression the integrator gets and tries to solve:
(%i7) :lisp (setq $form
'((MTIMES SIMP)
3* ((MEXPT SIMP)
((MPLUS SIMP) 1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))) -1)
1* ((MEXPT SIMP) $X -1)
2* ((MEXPT SIMP) $%E
((MTIMES SIMP) -1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))))))
If the routine intform finds a mexpt expression with a complicated
exponent it tries the following code:
(let* (($%e_to_numlog t)
(nexp (resimplify exp))) ; resimplify the whole integrand
(cond ((alike1 exp nexp) nil) ; no change, return nil
(t
(when *debug-integrate*
(format t "~& in INTFORM: NEW EXP is ~A~%" nexp))
(intform (setq exp nexp)))))))) ; try again with new exp
This code is involved doing the integral which fails. Now what happens
is that we get for every resimplify a new expression. The condition
(alike1 exp nexp) is never true and the integrator loops endlessly.
I have tried it directly at the Maxima prompt. Here are the results:
This is a first resimplify for the expression from above. Look at the
numbers:
(%i7) :lisp (resimplify $form)
((MTIMES SIMP)
2* ((MEXPT SIMP) $%E
((MTIMES SIMP) -1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))))
3* ((MEXPT SIMP)
((MPLUS SIMP) 1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))) -1)
1* ((MEXPT SIMP) $X -1))
We resimplify two times:
(%i7) :lisp (resimplify (resimplify $form))
((MTIMES SIMP)
1* ((MEXPT SIMP) $X -1)
2* ((MEXPT SIMP) $%E
((MTIMES SIMP) -1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))))
3* ((MEXPT SIMP)
((MPLUS SIMP) 1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))) -1))
And again:
(%i7) :lisp (resimplify (resimplify (resimplify $form)))
((MTIMES SIMP)
3* ((MEXPT SIMP)
((MPLUS SIMP) 1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1))) -1)
1* ((MEXPT SIMP) $X -1)
2* ((MEXPT SIMP) $%E
((MTIMES SIMP) -1 ((%LAMBERT_W SIMP) ((MEXPT SIMP) $X -1)))))
It seems to me that we have got by accident an expression which does not
simplify as expected and change the form every time we do a further
resimplify.
At the moment I have no idea how to continue.
Dieter Kaiser
_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima