#2732 wrong answer for similar to gaussian integral

None
closed
Dan Gildea
None
5
2014-06-04
2014-05-17
dan hayes
No

wxMaxima version: 13.4.0
Maxima version: 5.31.1
Maxima build date: 2013-09-24 09:49:12
Host type: i686-pc-mingw32
Lisp implementation type: GNU Common Lisp (GCL)
Lisp implementation version: GCL 2.6.8

(assume(a>0,b>0), t1:integrate(exp(-a^2*x^2-b^2/x^2),x,0,inf), t2:integrate(exp(-a^2*x^2-b^2/x^2),x,-inf,0)
, t3:integrate(exp(-a^2*x^2-b^2/x^2),x,-inf,inf),ldisp([t1,t2,t3]));

By symmetry of x and -x the first two results, t1 and t2 , should be same which they correctly are. But then obviously the 3rd result, t3, should just be twice t1 (or t2) but instead it adds a term with the the positive exponent in the exponential to t1 which is incorrect. Also maxima should be able to do closed form integrals of the form integrate(x^(2*i)*exp(-a*x^2-b/x^2),x,0,inf) for any integer i by differentiating i times the closed form answer of

integrate(exp(-a*x^2-b/x^2),x,0,inf)

with respect to a for i>0 or b for i<0 but it does not do so.

Discussion

  • For reference, here's the missing transcript (with current git master):

    (%i1) assume(a>0, b>0);
    (%o1)                           [a > 0, b > 0]
    (%i2) t1:integrate(exp(-a^2*x^2-b^2/x^2),x,0,inf);
                                              - 2 a b
                                  sqrt(%pi) %e
    (%o2)                         -------------------
                                          2 a
    (%i3) t2:integrate(exp(-a^2*x^2-b^2/x^2),x,-inf,0);
                                              - 2 a b
                                  sqrt(%pi) %e
    (%o3)                         -------------------
                                          2 a
    (%i4) t3:integrate(exp(-a^2*x^2-b^2/x^2),x,-inf,inf);
                                    2 a b               - 2 a b
                        sqrt(%pi) %e        sqrt(%pi) %e
    (%o4)               ----------------- + -------------------
                               2 a                  2 a
    (%i5) print([t1,t2,t3]);
                 - 2 a b              - 2 a b
     sqrt(%pi) %e         sqrt(%pi) %e
    [-------------------, -------------------, 
             2 a                  2 a
                                                      2 a b               - 2 a b
                                          sqrt(%pi) %e        sqrt(%pi) %e
                                          ----------------- + -------------------] 
                                                 2 a                  2 a
                       - 2 a b              - 2 a b
           sqrt(%pi) %e         sqrt(%pi) %e
    (%o5) [-------------------, -------------------, 
                   2 a                  2 a
                                                       2 a b               - 2 a b
                                           sqrt(%pi) %e        sqrt(%pi) %e
                                           ----------------- + -------------------]
                                                  2 a                  2 a
    
     
    • dan hayes
      dan hayes
      2014-05-24

      Note the commercial version Macsyma gives the correct answer of
      sqrt(%pi)exp(-2ab)/a for integrate(exp(-a^2x^2-b^2/x^2),x,-inf,inf);

      I see that the source code for that now defunct project is made publicly available so perhaps those knowledgeable could use or find out how Macsyma did it.

       
  • Dan Gildea
    Dan Gildea
    2014-06-04

    Fixed by [96309ddcd0dd398aade5f01eb1006512c23c2cf4]

    src/defint.lisp:
    o discontinuities-denom:
      recursively traverse exp in order to find discontinuities such as
      erfc(1/x-x) at x=0
    
    tests/rtestint.mac:
    o add integrate(exp(-x^2-1/x^2),x,-inf,inf);
    
     

    Related

    Commit: [96309d]

  • Dan Gildea
    Dan Gildea
    2014-06-04

    • status: open --> closed
    • assigned_to: Dan Gildea