NIntegrate giving message NIntegrate::slwcon:

Question

I got this interesting answer from Mathematica when trying to integrate my function numerically:

f[x_] := Sqrt[17*x^2 + x^4]
NIntegrate[f[x], {x, -1, 2}]

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

10.8053

But this is not an issue when integrating within other ranges. Why did I got this output?

Answer 1

I think the moral is that when the integrand has a singularity (its derivative is discontinuous at x = 0 in this case) and NIntegrate does not detect, the user has to tell it explicitly. That the singularity occurs when 0 is exactly 1/3 the way from an endpoint is interesting (in fact it can be off by a little), and I have a guess about how that slows convergence.

One standard way to handle a singularity is to add it to the domain specification like this:

NIntegrate[f[x], {x, -1, 0, 2}]

When this is done, the integral evaluates without a warning message to a precision of about 15 digits.

We can glimpse part of the problem by comparing Sqrt[17 x^2 + x^4] to an equivalent expression Abs[x] Sqrt[17 + x^2]. This last expression is not treated equivalently by NIntegrate. NIntegrate recognizes Abs as a piecewise function and divides the interval at 0, much as adding 0 to the domain above does. One can also turn off the symbolic processing by NIntegrate, in which case it evaluates the integral exactly as it does the OP's function.

To see, evaluate these and note the identical sampling and results:

{val, {pts}} = 
  Reap @ NIntegrate[Sqrt[17 x^2 + x^4], {x, -1, 0, 2}, EvaluationMonitor :> Sow[x]];
{val2, {pts2}} = 
  Reap @ NIntegrate[Abs[x] Sqrt[17 + x^2], {x, -1, 2}, EvaluationMonitor :> Sow[x]];
val == val2
pts == pts2
Length[pts]
(val - exact)/exact
(*
  True
  True
  22
  1.31517*10^-15
*)

{val, {pts}} = 
 Reap @ NIntegrate[Sqrt[17 x^2 + x^4], {x, -1, 2}, EvaluationMonitor :> Sow[x]];
{val2, {pts2}} = 
 Reap @ NIntegrate[Abs[x] Sqrt[17 + x^2], {x, -1, 2}, 
   EvaluationMonitor :> Sow[x], Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}];
val == val2
pts == pts2
Length[pts]
(val - exact)/exact

NIntegrate::slwcon: ...

NIntegrate::slwcon: ...

(*
  True
  True
  187
  -1.16022*10^-8
*)

One might note that the results are exactly the same in both cases, and especially that when singularity handling is turned on (in the first case), the result is found much more precisely with far fewer function evaluations.

Aside

The recursive subdivision carried out by NIntegrate divides the interval in half and resamples each half. If the error estimate in one half is too large, it will bisect it again. This continues until the error is acceptable or the MaxRecursion limit is reached. A number that is 1/3 the length of the interval from an endpoint is difficult to reach by this bisection process as it always stays 1/3 the length of the new subinterval from one of its endpoints. That cannot be the whole story, because it seems rare for a function to have this problem. Somehow the combination of the OP's particular function with the singularity being 1/3 the distance leads to the slow convergence of the numerical integration.

Here is the sampling on the OP integral:

Graphics[{
  Thin, Line[{{-1, 0}, {2, 0}}], PointSize[Medium],
  MapIndexed[{Hue[First[#2]/9], 
     Point[Function[x, {x, First[#2]}] /@ #1]} &, Split[pts, Less]]
  },
 AspectRatio -> 1/4, Frame -> True]

Update - Analysis of error

The tutorial subsection Example Implementation of a Global Adaptive Strategy shows how to analyze the error estimate. The code below is adapted from the tutorial. The default integration rule is "GaussKronrod" with 5 Gauss points and machine precision. IRuleEstimate below returns the approximation of the integral and an estimate of the error.

n = 5; precision = MachinePrecision;
{absc, weights, errweights} = NIntegrate`GaussKronrodRuleData[n, precision];

IRuleEstimate[f_, {a_, b_}] :=
 Module[{integral, error},
  {integral, 
    error} = (b - a) Total@
     MapThread[{f[#1] #2, f[#1] #3} &, {Rescale[absc, {0, 1}, {a, b}],
        weights, errweights}];
  {integral, Abs[error]}
  ]

We can see by comparing the estimate of the integral over the whole interval {-1, 2} and the interval broken up at 0, that there is tremendous difference in the value and error estimate. The recursive refinement used in the default "GlobalAdaptive" strategy will have a long way to go.

IRuleEstimate[f, {-1, 2}]
(* {10.7971, 0.224107} *)

IRuleEstimate[f, {-1, 0}] + IRuleEstimate[f, {0, 2}]
(* {10.8053, 6.52955*10^-10} *)

The "GlobalAdaptive" strategy subdivides the interval with the greatest error estimate. Here is code that sets up the first step step[0] and iterates the subdivision one step at a time.

setupGlobalAdaptive[f_, {a_, b_}] := Module[{t, integral, error},
   n = 5; precision = MachinePrecision;
   {absc, weights, errweights} = 
    NIntegrate`GaussKronrodRuleData[n, precision];

   {integral, error} = IRuleEstimate[f, {a, b}];
   {{{a, b}, integral, error}}
   ];

iterateGlobalAdaptive[f_, regions_] :=
  Module[{t, integral, error, r1, r2, a, b, c},

   (* splitting of the region with the largest error *)
   {a, b} = regions[[1, 1]]; c = (a + b)/2;

   (* integration of the left region *)
   {integral, error} = IRuleEstimate[f, {a, c}];
   r1 = {{a, c}, integral, error};

   (* integration of the right region *)
   {integral, error} = IRuleEstimate[f, {c, b}];
   r2 = {{c, b}, integral, error};

   (* sort the regions: the largest error one is the first *)
   Sort[
    Join[{r1, r2}, Rest[regions]],
    #1[[3]] > #2[[3]] &]
   ];

After eight further subdivisions of the integral over {-1, 2}, we arrive at the result of NIntegrate. One can see from the plot above or from the steps step[i], i = 0, 1, ... 8, that the subinterval with the greatest error contains the singular point 0.

step[0] = setupGlobalAdaptive[f, {-1, 2}];
Table[step[i] = iterateGlobalAdaptive[f, step[i - 1]]; 
 Total[Rest /@ step[i]], {i, 8}]
val - %[[-1, 1]]
(*
  {{10.8033, 0.0561962}, {10.8048, 0.0140596}, {10.8052, 0.00351555},
   {10.8053, 0.00087893}, {10.8053, 0.000219735}, {10.8053, 0.0000549339},
   {10.8053, 0.0000137335}, {10.8053, 3.43338*10^-6}}

  0.
*)

The Gauss-Kronrod rule performs well on smooth functions and much is known about its error for such functions. Less is known about non-smooth functions, but in practice, it performs well on many of them. For some reason the Gauss-Kronrod rule does not perform well on this function when 0 lies 1/3 from an endpoint in the interval of integration.

For the purposes of determining "slow convergence," the error is measured as a proportion of the value of the integral. While the relative error is persistently high for some domains of the form {-a, 2a} (as noticed on by Stephen Lutrell), it eventually stops complaining when the interval is large enough, for instance, when the interval is {-100, 200}.

There a couple of interesting things to note. One I have mentioned is the importance of singularities and helping NIntegrate to find them when it cannot do so automatically. The other is related. It is the importance of "SymbolicProcessing" in the case of Abs[x]. It helps NIntegrate above, although turning it off it is often suggested as a way to improve NIntegrate. All in all, I would say the mystery of this problem lies more in the mathematics than in Mathematica. A detailed analysis of the function and rule would probably reveal why the error is so great, but I leave that to anyone who has the leisure for it.

Answer 2

This is not an answer but a extended comment.

No more than Daniel Lichtblau, I don't know the why this is happening, but I have found four rational numbers in the range {-1, 3}, that produce the indicated message and for which small perturbations suppress the message.

With[{ϵ = 0}, NIntegrate[Sqrt[17 x^2 + x^4], {x, -1, #} + ϵ] & /@ {1/5, 1/2, 7/5, 2}]

produces NIntegrate::slwcon: for all four values.

With[{ϵ = 1*^-8}, NIntegrate[Sqrt[17 x^2 + x^4], {x, -1, #} + ϵ] & /@ {1/5, 1/2, 7/5, 2}]

produces NIntegrate::slwcon: for the first two values.

With[{ϵ = 1*^-7}, NIntegrate[Sqrt[17 x^2 + x^4], {x, -1, # + ϵ}] & /@ {1/5, 1/2, 7/5, 2}]

does not produce NIntegrate::slwcon: for any value. The behavior of Sqrt[17 x^2 + x^4 doesn't seem to be special in any way at these four points, so I am left puzzled.

I think this must be considered a bug.