Contradiction between Analytic and Numerical Integration

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).

• Whoops - I accidentally put two $dx$'s in my integral. I changed it now. Nov 15, 2016 at 22:01
• I suspect NIntegrate has trouble figuring out where the integrand is concentrated. But I confess that's just a guess. Nov 15, 2016 at 22:22

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  *)


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  *)


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

• Thanks. One question though - what does Mathematica do different in this case as compared to when a single Nintegrate[] was used? Nov 15, 2016 at 22:18
• not sure, but there are many situations where that works. Nov 15, 2016 at 22:21

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]
`