Contradiction between Analytic and Numerical Integration

Question

I have the function

f[x_, p_, b_] =  Exp[-p^2 - x^2]*Exp[4 b*p^2/(1 + 4 b)]/Pi/Sqrt[1 + 4 b]

and I want to find

\begin{align} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{p^2}{2}f(x,p,cix)dxdp \end{align} Integrating analytically, with $b=cix$, where $c\in\Re$ yields

In: Integrate[ Integrate[f[x, p, c*I*x]*(p^2/2), {p, -Infinity, Infinity},  Assumptions -> {c \[Element] Reals && x \[Element] Reals}], {x, -Infinity, Infinity}]
Out: 1/4

In contrast, using numerical integration, with $c=1$

In: NIntegrate[f[x, p, 1*I*x]*(p^2/2), {x, -Infinity, Infinity}, {p, -Infinity, Infinity}, MaxRecursion -> 20]
Out:0.266064 - 3.81639*10^-17 I

Or with $c=2$,

In: NIntegrate[f[x, p, 2*I*x]*(p^2/2), {x, -Infinity, Infinity}, {p, -Infinity, Infinity}, MaxRecursion -> 20]
Out:-0.0133697 + 6.00648*10^-17 I

Why are the numerical solutions so different than the analytic solution? I tried increasing the working precision to $20$, but I got the same answer (though with more decimal places).

Answer 1

I think the Levin Rule is messing up somehow (the default choice of Method here). I recall seeing such a mistake somewhere before. Here's a workaround:

NIntegrate[
 f[x, p, 1*I*x]*(p^2/2), {p, -Infinity, Infinity}, {x, -Infinity, Infinity}, 
 Method -> {"CartesianRule", Method -> {"GaussKronrodRule", "Points" -> 11}}]

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

(*  0.25 + 0. I  *)

---

Addendum

A way to get the Levin Rule to work is to make the substitution p^2 -> p and use the even symmetry of f in p to change the p integral to the domain {p, 0, Infinity}. This transforms the oscillatory part of the integrand into Exp[I * <linear>], which the Levin Rule seems to handle properly.

ff[x_, p_, b_] = Exp[-p - x^2]*Exp[4 b*p/(1 + 4 b)]/Pi/Sqrt[1 + 4 b];
NIntegrate[
 ff[x, p, 1*I*x]*(Sqrt[p]/2),
 {p, 0, Infinity}, {x, -Infinity, Infinity}]
(*  0.25 - 5.07359*10^-13 I  *)

Answer 2

If we nest the NIntegrate it gets it right..

 c=2;
 g[x_?NumericQ] := 
        NIntegrate[f[x, p, c*I*x]*(p^2/2), {p, -Infinity, Infinity}]
 NIntegrate[g[x], {x, -Infinity, Infinity}]

0.25 + 0. I

however.. reversing the order (which shouldn't matter) fails:

h[p_?NumericQ] := 
 NIntegrate[f[x, p, c*I*x]*(p^2/2), {x, -Infinity, Infinity}]
 NIntegrate[h[p], {p, -Infinity, Infinity}]

(convergence warning)

0.665207 - 0.0424908 I

Answer 3

The problem is quite easy to solve: The error arises, because NIntegrate does some symbolic pre-Processing of the integrand. In this case here, this preprocessing introduced a branch cut and the following numeric procession part got wrong results.You can see, that the funtion h has a jump at p= +-4.3. It can be avoided by setting SymbolicPreProcessing to False. This works both on the multivariate integral and the nested integrals.

   In[164]:= NIntegrate[
  f[x, p, 1*I*x]*(p^2/2), {x, -Infinity, Infinity}, {p, -Infinity, 
  Infinity}, Method -> {Automatic, "SymbolicProcessing" -> False}]

 During evaluation of In[164]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

  Out[164]= 0.25\[VeryThinSpace]- 2.49312*10^-16 I

 In[165]:= 
  h[p_?NumericQ] := 
 NIntegrate[f[x, p, 1*I*x]*(p^2/2), {x, -Infinity, Infinity}]
    NIntegrate[h[p], {p, -Infinity, Infinity}]

  During evaluation of In[165]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

  During evaluation of In[165]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in p near {p} = {-5.09582}. NIntegrate obtained 0.262508\[VeryThinSpace]-0.0736375 I
  and 0.003683811726300084` for the integral and error estimates. >>

    Out[166]= 0.262508\[VeryThinSpace]- 0.0736375 I

 In[167]:= 
  h2[p_?NumericQ] := 
  NIntegrate[f[x, p, 1*I*x]*(p^2/2), {x, -Infinity, Infinity}, 
   Method -> {Automatic, "SymbolicProcessing" -> False}]
  NIntegrate[h2[p], {p, -Infinity, Infinity}]

    Out[168]= 0.25\[VeryThinSpace]- 1.29797*10^-11 I

  In[169]:= Plot[Re@h[p], {p, -10, 10}, PlotRange -> All]

    In[170]:= Plot[Re@h2[p], {p, -10, 10}, PlotRange -> All]

  In[171]:= Plot[Im@h[p], {p, -10, 10}, PlotRange -> All]

  In[172]:= Plot[Im@h2[p], {p, -10, 10}, PlotRange -> All]