NIntegrate: Random occurence of “Catastrophic loss of precision” error

I have the following integral:

kernel[x_, xp_, phi_, d_] := 2 Sqrt[x^2 + xp^2 - 2 x xp Cos[phi]] - 2 Sqrt[d^2 + x^2 + xp^2 - 2 x xp Cos[phi]] - 1/2 d (Log[x^2 + xp^2 - 2 x xp Cos[phi]] + Log[-d + Sqrt[d^2 + x^2 + xp^2 - 2 x xp Cos[phi]]] - 3 Log[d + Sqrt[d^2 + x^2 + xp^2 - 2 x xp Cos[phi]]])
divM[x_,x0_] := (-1+(4 Sinh[x]^2)/(Cosh[2 x]+Cosh[2 x0])) (Sech[x-x0]+Sech[x + x0])+4/x Cosh[x0] Sinh[x]/(Cosh[2 x]+Cosh[2 x0])
Iv[x0_,d_] := NIntegrate[x xp divM[x,x0] divM[xp,x0] kernel[x,xp,phi,d], {x,0,Infinity}, {xp,0,Infinity}, {phi,0,Pi}, PrecisionGoal -> 3, Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 1000000}, MaxRecursion -> 20]

Or, in a human readable form:

$k(\varrho,\varrho^\prime,\phi,\delta)=-2 \sqrt{\delta ^2+\varrho ^2-2 \varrho \varrho^\prime \cos (\phi)+\varrho^{\prime 2}}+2 \sqrt{\varrho ^2-2 \varrho \varrho^\prime \cos (\phi)+\varrho^{\prime 2}}-\frac{1}{2} \delta \Big[\log \left(\sqrt{\delta ^2+\varrho ^2-2 \varrho \varrho^\prime \cos (\phi)+\varrho^{\prime 2}}-\delta \right)-3 \log \left(\sqrt{\delta ^2+\varrho ^2-2 \varrho \varrho^\prime \cos (\phi)+\varrho^{\prime 2}}+\delta \right)+\log \left(\varrho ^2-2 \varrho \varrho^\prime \cos (\phi)+\varrho^{\prime 2}\right)\Big]$

$\partial_\varrho(\varrho m_r(\varrho))= \frac{4 \cosh (\rho_0) \sinh (\varrho )}{\cosh (2 \rho_0)+\cosh (2 \varrho )}+\varrho \left[\frac{1}{\cosh(\varrho-\rho_0)}+\frac{1}{\cosh(\varrho+\rho_0)}\right] \left(\frac{4 \sinh ^2(\varrho )}{\cosh (2 \rho_0)+\cosh (2 \varrho )}-1\right)$

$I_v = \int_0^\infty d \varrho \int_0^\infty d \varrho^\prime \int_0^{\pi} d \phi \, \partial_\varrho(\varrho m_r(\varrho)) \partial_{\varrho^\prime}(\varrho^\prime m_r(\varrho^\prime)) k(\varrho,\varrho^\prime,\phi,\phi^\prime,\delta)$

Now my problem is that the function Iv fails to evaluate at random values of x0, for instance for x0=50. The error message is

Catastrophic loss of precision in the global error estimate due to insufficient WorkingPrecision or divergent integral.

Increasing the working precision has no effect. Also, there is no problem evaluating the integral for x0=40 or x0=60.

Any ideas what might cause this issue?

Also, for large values of x0>200, the integral fails when I use integer values for x0 (for instance Iv[1000,1]) and calculates without complaint, but to a wrong value, for floating point values (i.e., Iv[1000.,1]). Here, I am particularly worried about the wrong result that comes without a warning.

Edit:

As suggested in the comments, I tried evaluating the integral without the extra options, i.e.,

It[x0_,d_] := NIntegrate[x xp divM[x,x0] divM[xp,x0] kernel[x,xp,phi,d], {x,0,Infinity}, {xp,0,Infinity}, {phi,0,Pi}]

For x0=50, both It[50.,1] and Iv[50.,1] fail to evaluate. At x0=40, both evaluate, but yield very different results:

Iv[40.,1] = 16.86

It[40.,1] = 89.60

where the latter also says

The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 89.5960211577929 and 0.008059953558893987 for the integral and error estimates.

• 1. You are requesting a specific integration method, precision goal, etc. Just for kicks and giggles, have you tried to remove all your options and rerun the integration that way? Maybe the problem is in the interaction of two options. 2. Have you tried evaluating the functions in the Iv integrand on their own, to see what their values are when you get that loss of precision? It would help to perhaps pinpoint the origin of the problem. – MarcoB Jun 1 '16 at 14:29
• For 1, see my edits. For 2: The difference between x0=40 and x0=50 is basically that divM is shifted by 10 units along the x-axis (in a very good approximation, divM[x,50] = divM[x-10,40]. Nothing special happens at x0=50. – Felix Jun 1 '16 at 14:50

You are using smaller precision goal when computing with Iv, that is why you get a different result. If you increase the MinRecursion value you get the same result as with It:

In:= Clear[Iv]
Iv[x0_, d_] :=
NIntegrate[
x xp divM[x, x0] divM[xp, x0] kernel[x, xp, phi, d], {x, 0,
Infinity}, {xp, 0, Infinity}, {phi, 0, Pi}, PrecisionGoal -> 3,
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 1000000},
MinRecursion -> 4, MaxRecursion -> 20]

In:= Iv[40., 1]

Out= 89.8663

And (only) the message:

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.

• You are right. I was assuming that a precision goal of 3 kind of guarantees that the result is accurate to 3 digits. I came to realize that this is not the case and now always double check if the result of NIntegrate (or basically any other computation) is reasonable. – Felix Jan 27 '17 at 15:41