20

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.

Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
DumpsterDoofus
  • 11,857
  • 1
  • 29
  • 49
  • Try replacing Sqrt[x] with x^0.5 – bill s Nov 03 '13 at 20:27
  • 2
    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 – Nasser Nov 03 '13 at 21:00
  • The integral works correctly in V7. It works without Assumptions in V8.0.4, but not with the option. – Michael E2 Nov 03 '13 at 21:23
  • @Nasser: Then aren't both Mathematica and Maple giving wrong answers? The correct answer is 0.907, so irregardless of whether there is an $i$ or not, they are both incorrect. – DumpsterDoofus Nov 03 '13 at 21:26
  • Sure. So both Maple and Mathematica are wrong. For some reason, now I feel a little better :) – Nasser Nov 03 '13 at 21:32
  • @Nasser The $i$ is most easily explained by a different choice of branches for the fourth root in (a^2)^(1/4), which we see in your Maple output. – Michael E2 Nov 03 '13 at 21:33
  • @MichaelE2 good points about branch cuts. But I have thought the math/computer expert folks have decided on these things long time ago (PrincipalValue?) else one CAS will give one result, depending on its choice, and another will give different result for same problem because it decided to go another path (as it seems in this case). Any way, will have to learn more about these things myself one day. Good observations, thank you. – Nasser Nov 03 '13 at 22:12
  • 3
    @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^(2ipi) 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). – asmeurer Nov 04 '13 at 01:42
  • 6
    It's a bug introduced in version 8. We're looking into it. – Daniel Lichtblau Nov 04 '13 at 22:59

1 Answers1

16

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 *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thanks for replying! That definitely works. However, I still don't understand why the original code is giving contradictory and incorrect answers. – DumpsterDoofus Nov 03 '13 at 20:34
  • @DumpsterDoofus Me neither, as yet. It certainly looks like something to do with branch cuts, but I haven't figured it out. :/ – Michael E2 Nov 03 '13 at 20:40
  • 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. – DumpsterDoofus Nov 03 '13 at 21:28
  • This the same answer you get using x^0.5 instead of Sqrt[x], so this is the likely culprit. – bill s Nov 03 '13 at 21:36