Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
Try replacing Sqrt[x] with x^0.5 –  bill s Nov 3 '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 3 '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 3 '13 at 21:23
2  
@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). –  asmeurer Nov 4 '13 at 1:42
5  
It's a bug introduced in version 8. We're looking into it. –  Daniel Lichtblau Nov 4 '13 at 22:59

1 Answer 1

up vote 13 down vote accepted

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 *)
share|improve this answer
    
Thanks for replying! That definitely works. However, I still don't understand why the original code is giving contradictory and incorrect answers. –  DumpsterDoofus Nov 3 '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 3 '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 3 '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 3 '13 at 21:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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