Numerical integral of finite valued function over finite integral fails with NIntegrate::inumri

I have the following (simplified) integral:

f[x_, x0_] := Sinh[x]/(Cosh[2 x] + Cosh[2 x0])
g[x_, xp_] := EllipticK[(4 x xp)/(x + xp)^2]/(x + xp)
NIntegrate[
f[x, 101.] f[xp, 101.] g[x, xp], {x, 51., 101.}, {xp, 50., x-10.^-6},
PrecisionGoal -> 6, MaxRecursion -> 100]


For me (with Mathematica v. 10.3), the integral fails with

NIntegrate::inumri: "The integrand (...) has evaluated to Overflow, Indeterminate, or Infinity for all sampling points the region with boundaries {{0.,9.35762*10^-14},{0,1}}"

Why is that? Both f and g are finite and well defined for all values (in particular in the integration range).

Also, the "boundaries" given in the error message do not make sense to me. How do I read this? None of these values are in the interval [51, 101].

When I omit the MaxRecursion parameter, the integral gives a warning instead of an error saying that it could not get a correct value with the given number of recursions (wrong number and no number seem equally bad to me).

Edit:

It also fails when defining f and g to be numeric:

f[x_?NumericQ, x0_?NumericQ] := Sinh[x]/(Cosh[2 x] + Cosh[2 x0])
g[x_?NumericQ, xp_?NumericQ] :=
EllipticK[(4 x xp)/(x + xp)^2]/(x + xp)
NIntegrate[
f[x, 101.] f[xp, 101.] g[x, xp], {x, 51., 101.}, {xp, 50.,
x - 10.^-6},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000},
MaxRecursion -> 100]


The MaxErrorIncreases parameter was required to avoid warnings about inaccurate results.

In my real application, f has an additional factor Cosh[x0] that I omitted here because it is a constant factor in the integral. Just now, it happens that f evaluates to very small numbers, but again, this should only be a constant factor. The same result is obtained with WorkingPrecision -> 20.

Edit2:

Since there was some concern about the small values of the function f, here another function that gives the same problem:

f2[x_?NumericQ, x0_?NumericQ] := Sech[x - x0] (1 - x Tanh[x - x0])
NIntegrate[
f2[x, 100.] f2[xp, 100.] g[x, xp], {x, 51., 101.}, {xp, 50.,
x - 10.^-6},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000},
MaxRecursion -> 100]


With the previous value of x0=101, this new integral gives a value without any error messages, but with x0=100, it fails as before. It seems to be quite random when the integral can be calculated and when not.

• Note that g[51, 51], which happens to be on one of the boundaries of your region of integration, evaluates to ComplexInfinity, so perhaps your functions are not as finite-valued as you thought. – MarcoB Jun 9 '16 at 16:32
• That's true, well spotted. However, I get the same error when excluding the point x=xp from the integration (i.e., by stopping at xp=x-1e-6, see edits to the question. – Felix Jun 9 '16 at 18:44
• You may want to prevent symbolic evaluation of f and g using NumericQ; that's typically safer for numerical evaluations: f[x_?NumericQ, x0_?NumericQ] := ... and similarly with g. See mathematica.stackexchange.com/a/26037/27951. I'm not sure that this will fix the problem though. The f and g functions evaluate to very very small numbers over your domain. Instead of a PrecisionGoal, try using only exact numbers in your definitions, and setting WorkingPrecision -> \$MachinePrecision (or higher) to enable arbitrary-precision arithmetic and error tracking. – MarcoB Jun 9 '16 at 20:10
• Thanks for this tip, but unfortunately, the problem persists (see edit). – Felix Jun 9 '16 at 20:53
• Have you tried plotting your integrand with Plot3D[], just to check that what you're integrating is not nasty-looking? – J. M.'s torpor Jun 9 '16 at 21:08

Let us have a single integrand function first:

F[x_?NumericQ, xp_?NumericQ] := (
EllipticK[(4 x xp)/(x + xp)^2] Sinh[x] Sinh[
xp])/((x + xp) (Cosh[202] + Cosh[2 x]) (Cosh[202] + Cosh[2 xp]));


obtained by expanding f[x, 101.] f[xp, 101.] g[x, xp].

For me (with Mathematica v. 10.3), the integral fails with

NIntegrate::inumri: "The integrand (...) has evaluated to Overflow, Indeterminate, or Infinity for all sampling points the region with boundaries {{0.,9.35762*10^-14},{0,1}}"

Why is that? Both f and g are finite and well defined for all values (in particular in the integration range).

This is because of the singularity handler "DuffyCoordinates". With no singularity handler or "IMT" the integration process finishes only issuing "NIntegrate::slwcon".

NIntegrate[F[x, xp], {x, 51., 101.}, {xp, 50., x - 10.^-6},
Method -> {"GlobalAdaptive", "SingularityHandler" -> None},
PrecisionGoal -> 6, MaxRecursion -> 100]

(* 2.06466*10^-90 *)


If larger precision is used again only the message "NIntegrate::slwcon" is issued:

NIntegrate[F[x, xp], {x, 51, 101}, {xp, 50, x - 10^-6},
WorkingPrecision -> 60, PrecisionGoal -> 6, MaxRecursion -> 100]

(*2.06466493230756122140054681757475975443145922411426227150693*10^-90 *)


Sampling points

Compare the following plots of integration sampling points. The first is with the automatic singularity handling, the second with IMT, the third without any:

Needs["IntegrationNIntegrateUtilities"]

Grid[{
Table[
(k = 0;
NIntegrate[F[x, xp], {x, 51., 101.}, {xp, 50., x - 10.^-6},
Method -> {"GlobalAdaptive", "SingularityHandler" -> sh},
EvaluationMonitor :> (k++), PrecisionGoal -> 6,
MaxRecursion -> 100];
gr = NIntegrateSamplingPoints[
NIntegrate[F[x, xp], {x, 51., 101.}, {xp, 50., x - 10.^-6},
Method -> {"GlobalAdaptive", "SingularityHandler" -> sh},
PrecisionGoal -> 6, MaxRecursion -> 100]];
Append[Append[gr, ImageSize -> Medium],
PlotLabel ->
Row[{"SingularityHandler: ", sh,
"\nnumber of sampling points: ", k}]]), {sh, {Automatic,
"IMT", None}}]}, Dividers -> All]


It seems that using IMT produces (seemingly) less points.

With larger working precision (60) much less sampling points are needed:

This can be explained with integrand's structure. Additionally, with "DuffyCoordinates" the integrand gets more complicated.

• That makes sense. Now it is also clear why the error message complains about a parameter range that is not within the integration range: Because of the transformation of variables by the singularity handler. – Felix Jun 10 '16 at 19:37
• Correct, I thought but forgot to mention that... – Anton Antonov Jun 10 '16 at 19:38
• How peculiar. Duffy is supposed to be the "cheaper" substitution, since it does not involve a transformation with transcendental equations unlike IMT. Also, with IMT, you are then forced to do a Cartesian product integration, which makes it unsuitable for Genz-Malik. But I guess the cost of evaluating the integrand itself trumps everything else. – J. M.'s torpor Jun 11 '16 at 4:23
• @J.M. This kind of seems to be the case -- see the updated plots in my answer. – Anton Antonov Jun 11 '16 at 14:35