0
$\begingroup$

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.

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Jun 9, 2016 at 16:32
  • $\begingroup$ 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. $\endgroup$
    – Felix
    Commented Jun 9, 2016 at 18:44
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Jun 9, 2016 at 20:10
  • $\begingroup$ Thanks for this tip, but unfortunately, the problem persists (see edit). $\endgroup$
    – Felix
    Commented Jun 9, 2016 at 20:53
  • $\begingroup$ Have you tried plotting your integrand with Plot3D[], just to check that what you're integrating is not nasty-looking? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 9, 2016 at 21:08

1 Answer 1

Reset to default
2
$\begingroup$

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["Integration`NIntegrateUtilities`"]

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]

enter image description here

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

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

enter image description here

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Felix
    Commented Jun 10, 2016 at 19:37
  • $\begingroup$ Correct, I thought but forgot to mention that... $\endgroup$
    – Anton Antonov
    Commented Jun 10, 2016 at 19:38
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 11, 2016 at 4:23
  • $\begingroup$ @J.M. This kind of seems to be the case -- see the updated plots in my answer. $\endgroup$
    – Anton Antonov
    Commented Jun 11, 2016 at 14:35

Your Answer

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

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