7
$\begingroup$

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

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

3 Answers 3

10
$\begingroup$

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  *)
$\endgroup$
5
$\begingroup$

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

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

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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