Integrate returns imaginary answer for smooth, real function

Question

Bug introduced in 7 or earlier and fixed in 14.0.0

---

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, 10.3 and 12.1 all give the same result.

Answer 1

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[238]= (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[239]= -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[270]= 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[266]= (E^-b^2 (-Log[-b] + Log[b]))/(4 b Sqrt[\[Pi]])
*)

(*
Out[267]= 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[273]= 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[328]= (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[330]= 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[292]= (
Gamma[1 - p/2] Hypergeometric1F1[(1 + p)/2, 3/2, -b^2])/(2 Sqrt[\[Pi]])
*)

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

(*
Out[294]= 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[307]= 1/4 (-2 \[Pi] (2 I + Erfi[1/2]) + HypergeometricPFQ[{1, 1}, {3/2, 2}, 1/4])
*)

(*
Out[308]= -0.693664 - 3.14159 I
*)

The true value is

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

(*
Out[313]= 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[310]= 1/2 (\[Pi] Erfi[1/2] + 1/2 HypergeometricPFQ[{1, 1}, {3/2, 2}, 1/4])
*)

(*
Out[311]= 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[320]= 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[320]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

(*
Out[321]= 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[323]= 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.

Answer 2

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)