Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.
Close
From: Raymond Toy <rtoy@us...>  20060327 22:57:31

Update of /cvsroot/maxima/maxima/src In directory sc8prcvs1.sourceforge.net:/tmp/cvsserv7940/src Modified Files: defint.lisp Log Message: o Add comments o Add LOGTRANSFORM to convert an integral of the form p(x)*f(exp(x)) for x from minf to inf to the equivalent p(log(y))*f(y) for y from 0 to inf. p(x) is a polynomial. o Adjust MTOINF for the case of p(x)*f(exp(x))  Say the integral diverges only if the limit of the integrand is not zero at both minf and inf. (Are there other cases?)  Add a call to LOGTRANSFORM in case RECTZTO%PI2 fails, which happens when f(z) is not a rational function. The transformation should be to a form that maxima knows how to handle.  If LOGTRANSFORM fails, just give up so we try other methods or return the integral itself. This should fix bug 1451351. Index: defint.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v retrieving revision 1.20 retrieving revision 1.21 diff u d r1.20 r1.21  defint.lisp 27 Mar 2006 17:34:52 0000 1.20 +++ defint.lisp 27 Mar 2006 22:57:26 0000 1.21 @@ 1293,13 +1293,25 @@ (involve grand '(%sinh %cosh %tanh))) (p*pin%ex n) ;setq's P* and PE*...Barf again. (setq ans (catch 'pin%ex (pin%ex d)))) + ;; We have an integral of the form p(x)*F(exp(x)), where + ;; p(x) is a polynomial. (cond ((null p*) + ;; No polynomial (return (dintexp grand var))) ((not (and (zerop1 (getlimit grand var '$inf)) (zerop1 (getlimit grand var '$minf)))) ;; These limits must exist for the integral to converge. (diverg)) ((setq ans (rectzto%pi2 (m*l p*) (m*l pe*) d)) + ;; This only handles the case when the F(z) is a + ;; rational function. + (return (m* (m// nc dc) ans))) + ((setq ans (logtransform (m*l p*) (m*l pe*) d)) + ;; If we get here, F(z) is not a rational function. + ;; We transform it using the substitution x=log(y) + ;; which gives us an integral of the form + ;; p(log(y))*F(y)/y, which maxima should be able to + ;; handle. (return (m* (m// nc dc) ans))) (t ;; Give up. We don't know how to handle this. @@ 2264,6 +2276,16 @@ (m (m+ n 1) (m+ m 1))) ((signp l m) (solvex eql cl nil nil)))) +;; Integrate(p(x)*f(exp(x))/g(exp(x)),x,minf,inf) by applying the +;; transformation y = exp(x) to get +;; integrate(p(log(y))*f(y)/g(y)/y,y,0,inf). This should be handled +;; by dintlog. +(defun logtransform (p pe d) + (let ((newp (subst (list '(%log) var) var p)) + (newpe (subst var 'z* (catch 'pin%ex (pin%ex pe)))) + (newd (subst var 'z* (catch 'pin%ex (pin%ex d))))) + (defint (div (div (mul newp newpe) newd) var) var 0 ul))) + ;; This implements Wang's algorithm in Chapter 5.2, pp. 98100. ;; ;; This is a very brief description of the algorithm. Basically, we @@ 2389,7 +2411,7 @@ ;;(declaretop (special e)) ;; Test to see if exp is of the form f(exp(x)), and if so, replace ;; exp(x) with 'z*. +;; exp(x) with 'z*. That is, basically return f(z*). (defun pin%ex (exp) (declare (special $exponentialize)) (pin%ex0 (cond ((notinvolve exp '(%sinh %cosh %tanh)) @@ 2399,7 +2421,10 @@ (setq exp ($expand exp))))))) (defun pin%ex0 (e)  (declare (special e)) ; Why does e need to be special? + ;; Does e really need to be special here? Seems to be ok without + ;; it; testsuite works. + #+nil + (declare (special e)) (cond ((not (among var e)) e) ((atom e) 