# Integrate returns imaginary answer for smooth, real function

Bug introduced in 7 or earlier and persisting through 12.0.0.0 or later

I'm trying to evaluate the integral:

$$\int_0^{\infty} \frac{1}{4 b \sqrt{\pi} r} e^{-(b-r)^2}(e^{4 b r} - 1) \mathrm{d}r$$

with $$b>0$$.

Integrate[(E^-(b + r)^2 (-1 +
E^(4 b r)))/(4 b Sqrt[π] r), {r, 0, ∞},
Assumptions -> b > 0]
(* -((I E^-b^2 Sqrt[π])/(4 b)) *)


Unfortunately, the function being integrated is purely real and smooth, so the answer should be real, but the one provided by Mathematica is purely imaginary.

-((I E^-b^2 Sqrt[π])/(4 b)) /. b -> 1/(π + E) // N
(* 0. - 2.52206 I *)


NIntegrate is able to correctly return real values.

NIntegrate[-(E^-(b + r)^2/(4 b Sqrt[π] r)) +
E^(4 b r - (b + r)^2)/(4 b Sqrt[π] r) /.
b -> 1/(π + E), {r, 0, ∞}]
(* 0.490405 *)


Unlike some other cases I have seen, Mathematica does not know the elementary antiderivative for this function, so I cannot tell whether some issue with discontinuities or branch cuts is causing this problem. What explains this strange result; is it a bug?

I've checked that versions 9.0, 10.2, and 10.3 all give the same result.

• note you can demonstrate the same issue with a specific value for b – george2079 Nov 13 '15 at 21:50
• Change of variables r -> 1/x gives the result as DawsonF[b]/(2b) – Simon Woods Nov 13 '15 at 22:09
• Reported as a bug. – Daniel Lichtblau Nov 13 '15 at 22:51
• It has given the same answer since version 7 at the latest. Version 5.2 gave a different answer, also wrong (it gives a complex result which is the sum of the correct, purely real result and the negative of the erroneous, purely imaginary one). – Oleksandr R. Nov 14 '15 at 23:12
• @Kevin Driscoll: Notice that the Latex expression differs from your code expression in the sign of b in the exponent. Now, whereas the code version gives your imaginary result, the Latex version gives a real result in terms of Erfi and HypergeometricPFQ which is, however, also wrong as it is negative for b>0 despite of the positive integrand. Hence we face a fairly buggy situation which extends even beyond your discovery. – Dr. Wolfgang Hintze Nov 18 '15 at 11:39

Although this question has been answered by a bug statement, I found it interesting to study some aspects in more detail.

The study could not explain the internal cause of the bug.

But, besides taking into account that the OP provides two different integrals the study shows that applying a simple unified regularization procedure leads to correct results.

The two integrals of the OP

The Latex version of the integral is

intLatex = Integrate[(Exp[-(b - r)^2] (-1 + Exp[4 b r]))/(4 b r Sqrt[\[Pi]]), {r, 0, \[Infinity]}, Assumptions -> b > 0]
% /. b -> 1/2.

(*
Out= (1/(8 b Sqrt[\[Pi]]))E^-b^2 (\[Pi] (Erfi[b] - Erfi[3 b]) -
2 b^2 (HypergeometricPFQ[{1, 1}, {3/2, 2}, b^2] -
9 HypergeometricPFQ[{1, 1}, {3/2, 2}, 9 b^2]))
*)

(*
Out= -0.1717
*)


It gives the wrong result of a negative value even though the integrand is positive.

For comparison, the numeric integration gives the correct value

intLatexN =
NIntegrate[(Exp[-(b - r)^2] (-1 + Exp[4 b r]))/(4 b r Sqrt[\[Pi]]) /.
b -> 1/2., {r, 0, \[Infinity]}]

(*
Out= 2.56822
*)


The version of the integral provided in the MMA code is (it has (r+b) instaed of (r-b) in the exponent

 intCode = Integrate[(E^-(b + r)^2 (-1 + E^(4 b r)))/(4 b Sqrt[\[Pi]] r), {r, 0, \[Infinity]}, Assumptions -> a < 2]
% /. b -> 1/2.

(*
Out= (E^-b^2 (-Log[-b] + Log[b]))/(4 b Sqrt[\[Pi]])
*)

(*
Out= 0. - 0.690194 I
*)


Again the numeric values for comparison:

intCodeN =
With[{b = 1/2},
NIntegrate[(E^-(b + r)^2 (-1 + E^(4 b r)))/(4 b Sqrt[\[Pi]] r), {r,
0, \[Infinity]}]]

(*
Out= 0.424436
*)


Hence the symbolic integration over the positive integrand leads to a purely imaginary value which has no appearent relation to the true value.

Modifications

Here we discuss two modifications, first we replace the inverse of $r$ by a real power of $r$, and secondly we provide an even simpler example of the buggy situation.

Denominator

If we modify the denominator of the integrands replacing 1/r by 1/r^p, where p > 0 is an auxiliary parameter, we get correct results in the limit p->1

intLatexp =
Integrate[(Exp[-(b - r)^2 ] (-1 + E^(4 b r)))/(4 b Sqrt[\[Pi]] r^p), {r,
0, \[Infinity]}, Assumptions -> p < 2]
Limit[%, p -> 1]
% /. b -> 1/2.

(*
Out= (1/(8 b Sqrt[\[Pi]]))E^-b^2 (Gamma[
1/2 - p/2] (-Hypergeometric1F1[1/2 - p/2, 1/2, b^2] +
Hypergeometric1F1[1/2 - p/2, 1/2, 9 b^2]) +
b p Gamma[-(p/2)] (Hypergeometric1F1[1 - p/2, 3/2, b^2] -
3 Hypergeometric1F1[1 - p/2, 3/2, 9 b^2]))
*)


$\frac{e^{-b^2} \left(-\text{Hypergeometric1F1}^{(1,0,0)}\left(0,\frac{1}{2},b^2\right)+\text{Hypergeometric1F1}^{(1,0,0)}\left(0,\frac{1}{2},9 b^2\right)+\pi (\text{erfi}(3 b)-\text{erfi}(b))\right)}{8 \sqrt{\pi } b}$

(*
Out= 2.56822
*)

intCodep =
Integrate[(Exp[-(b + r)^2 ] (-1 + E^(4 b r)))/(4 b Sqrt[\[Pi]] r^p), {r, 0, \[Infinity]}, Assumptions -> p < 2]
% /. p -> 1
% /. b -> 1/2.

(*
Out= (
Gamma[1 - p/2] Hypergeometric1F1[(1 + p)/2, 3/2, -b^2])/(2 Sqrt[\[Pi]])
*)

(*
Out= (E^-b^2 Sqrt[\[Pi]] Erfi[b])/(4 b)
*)

(*
Out= 0.424436
*)


Simplified problem

This integral has a similar problem: positive integrand but complex result with negative real part

intSim = Integrate[Exp[-r^2] (Exp[r] - 1)/r, {r, 0, \[Infinity]}]
% // N

(*
Out= 1/4 (-2 \[Pi] (2 I + Erfi[1/2]) + HypergeometricPFQ[{1, 1}, {3/2, 2}, 1/4])
*)

(*
Out= -0.693664 - 3.14159 I
*)


The true value is

intSimN = NIntegrate[Exp[-r^2] (Exp[r] - 1)/r, {r, 0, \[Infinity]}]

(*
Out= 1.23826
*)


This substitution removes the error

r^2 -> z, dr = 1/2 dz/Sqrt[z]

intSims = 1/2 Integrate[Exp[-z] (Exp[Sqrt[z]] - 1)/z, {z, 0, \[Infinity]}]
% // N

(*
Out= 1/2 (\[Pi] Erfi[1/2] + 1/2 HypergeometricPFQ[{1, 1}, {3/2, 2}, 1/4])
*)

(*
Out= 1.23826
*)


Again the real power limit gives the correct result:

intSimp = Integrate[Exp[-r^2] (Exp[r] - 1)/r^p, {r, 0, \[Infinity]},
Assumptions -> p < 2]
% /. p -> 1 (* does not work Limit required *)
Limit[%%, p -> 1]
% // N

(*
Out= 1/2 (Gamma[1/2 - p/2] (-1 + Hypergeometric1F1[1/2 - p/2, 1/2, 1/4]) +
Gamma[1 - p/2] Hypergeometric1F1[1 - p/2, 3/2, 1/4])
*)


During evaluation of In:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

(*
Out= Indeterminate
*)


$\frac{1}{2} \left(\text{Hypergeometric1F1}^{(1,0,0)}\left(0,\frac{1}{2},\frac{1}{4}\right)+\pi \text{erfi}\left(\frac{1}{2}\right)\right)$

(*
Out= 1.23826
*)


Conslusion

Although we could not clarify the reason for the bug we have shown that it appears similarly in even simpler related situations.

All integrals investigated here can be given a consitent meaning if understood as a limit p->1 of the integral over an integrand "regularized" by replacing 1/r by 1/r^p.

• Not so happy with somebody downvoting my contribution instead of at least taking care of my comment with respect to the inconsistency in the OP, which nevertheless received a rather high vote. But I have somehow given up to understand the voting behaviour in this forum. Finally, at least for me the answer was fun and useful too :-) – Dr. Wolfgang Hintze Nov 22 '15 at 20:31
• Stumbled upon this thread after starting to type a question about the same bug. Seems like MMA doesn't like anything more complicated than -x^2 in the exponent; check out my approach and +1 to you. – LLlAMnYP Nov 7 '17 at 14:22

I at first failed to find this question and started writing my own when I bumped into this bug. My MWE is even more trivial than Wolfgang's, so I'm posting this answer in hope that it might be an easier workaround for certain problems.

Integrate[Exp[-(x - x0)^2]/x, {x, -Infinity, Infinity},
PrincipalValue -> True,
Assumptions -> x0 > 0]

-I E^-x0^2 \[Pi]


The integrand is undoubtedly real-valued. The Cauchy principal value has every reason to exist and be real since taking

integrand[x_] := Exp[-(x - x0)^2]/x
Limit[integrand[x] + integrand[-x], x -> 0]

4 E^-x0^2 x0 (* finite real value *)


and

Plot[integrand[x] + integrand[-x] /. x0 -> 1, {x, 0, 4}] looks perfectly well-behaved. But I don't need to explain this here really.

This problem is trivially circumvented by an even easier "regularization".

Integrate[
Exp[-(x - x0)^2]/x /. x -> x + x0, {x, -Infinity, Infinity},
PrincipalValue -> True,
Assumptions -> x0 > 0]

2 Sqrt[π] DawsonF[x0]


By a simple coordinate transformation we got the correct answer.

The problem in OP, by the way, is almost the same one. Let's note that his integrand is simply

integrand1[r_] := (Exp[-(b-r)^2] - Exp[-(b+r)^2])/(4 b Sqrt[Pi] r)


so it's quite easy to simplify with

integrand2[r_] := Exp[-(b-r)^2]/(4 b Sqrt[Pi] r)
integrand1[r] == integrand2[r] + integrand2[-r] // FullSimplify

True


Therefore

Integrate[integrand1[r], {r, 0, Infinity}]


is equivalent to

Integrate[integrand2[r], {r, -Infinity, Infinity}]


which is the same as

Integrate[integrand2[r + b], {r, -∞, ∞}, PrincipalValue -> True, Assumptions -> b > 0]

DawsonF[b]/(2 b)