20
$\begingroup$

Bug introduced in 8.0 and fixed in 10.0


I attempted to calculate the following integral:

Integrate[ Sqrt[x] Exp[-(x - a)^2], {x, 0, ∞}, Assumptions -> a > 0]

1/(4 Sqrt[a])I E^(-(a^2/2)) π( a^2 BesselI[-(1/4), a^2/2] - (1 + a^2) BesselI[1/4, a^2/2]
                             + a^2 (BesselI[3/4, a^2/2] - BesselI[5/4, a^2/2]))

This is obviously incorrect, since applying the rule $a \rightarrow 0.3$ yields $-0.37\mathrm{i}$, whereas the correct result obtained from NIntegrate is $0.907$.

I have two questions:

  1. Does anyone know the correct answer? And perhaps more importantly,
  2. Does anyone have any idea why such a simple integral is resulting in a completely incorrect output? I find it disturbing that the program is capable of yielding incorrect results without warning to the user. In this case the answer was obviously wrong, but what about for more complicated expressions?

Note: removing the Assumptions option still yields the same incorrect result.

$\endgroup$
9
  • $\begingroup$ Try replacing Sqrt[x] with x^0.5 $\endgroup$
    – bill s
    Nov 3, 2013 at 20:27
  • 2
    $\begingroup$ Something is wrong. I tried this on Maple, with no assumptions. Mathematica and Maple gave the same numerical answer, except Mathematica made it imaginary (with minus sign), while Maple showed it as real (with positive sign). i.e. both gave 0.37 (in absolute terms). (again, no assumptions). So either Maple or Mathematica has a bug :) here is a screen shot: !Mathematica graphics $\endgroup$
    – Nasser
    Nov 3, 2013 at 21:00
  • $\begingroup$ The integral works correctly in V7. It works without Assumptions in V8.0.4, but not with the option. $\endgroup$
    – Michael E2
    Nov 3, 2013 at 21:23
  • 3
    $\begingroup$ @Nasser it's not that easy. You really have to keep track of what branch you're on throughout computations. For instance, if a subcomputation gives e^(2*i*pi) and replaces it immediately it with 1, then if a later computation takes the square root the answer will be wrong, because it should have been on the branch where sqrt(1) = -1 (essentially). $\endgroup$
    – asmeurer
    Nov 4, 2013 at 1:42
  • 6
    $\begingroup$ It's a bug introduced in version 8. We're looking into it. $\endgroup$ Nov 4, 2013 at 22:59

1 Answer 1

16
$\begingroup$

Here's the exact answer:

i1 = Integrate[x^n Exp[-(x - a)^2], {x, 0, Infinity}, 
   Assumptions -> n > 0] /. n -> 1/2

(*  1/2 E^-a^2 (Gamma[3/4] Hypergeometric1F1[3/4, 1/2, a^2] + 
      1/2 a Gamma[1/4] Hypergeometric1F1[5/4, 3/2, a^2])  *)

i1 /. a -> 0.3
(* 0.907605 *)
$\endgroup$
4
  • $\begingroup$ Thanks for replying! That definitely works. However, I still don't understand why the original code is giving contradictory and incorrect answers. $\endgroup$ Nov 3, 2013 at 20:34
  • $\begingroup$ @DumpsterDoofus Me neither, as yet. It certainly looks like something to do with branch cuts, but I haven't figured it out. :/ $\endgroup$
    – Michael E2
    Nov 3, 2013 at 20:40
  • $\begingroup$ Yeah, I figured it might have something to do with how it's handing branch cuts, but I don't know enough about how the MathKernel works to guess any further. $\endgroup$ Nov 3, 2013 at 21:28
  • $\begingroup$ This the same answer you get using x^0.5 instead of Sqrt[x], so this is the likely culprit. $\endgroup$
    – bill s
    Nov 3, 2013 at 21:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.