# How much can we trust NIntegrate when NIntegrate::eincr is shown?

I have been dealing with a high oscillation integration using NIntegrate. The warning of eincr shows that no matter how large (X) I use in the MaxErrorIncreases option, "global error of the strategy GlobalAdaptive has increased more than X times". Or rather, I cannot get rid of this warning, although the estimated error given by the warning is indeed decreasing.

Then I tried to use the Cuba package, which is used to evaluate high oscillation integrations in QCD. (http://www.feynarts.de/cuba/) The evaluation process is very slow and the estimated error after Maxpoints -> 10000 evaluation is still larger than the estimated error from NIntegrate.

Should I believe in NIntegrate more than Cuba in this case? Or both are just bad results?

Sorry I don't think it is a good idea to post the total code online. My supervisor will be pissed off. The project is doing a large 5 dimensional integration over two phase spaces that are Lorentz transformed two times. The integrand is that quantity which mix up two phase spaces' variables by an arbitrary Lorentz transformation two times, which contains hundreds of terms. I believe that the project cannot be simplified more. I just want to know the meaning of "global error of the strategy GlobalAdaptive has increased more than X times".

Below are the warnings. I believe they are commonly seen.

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.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 100000 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 -0.00027648 and 0.00003375392679593996 for the integral and error estimates.

• Not much to say without your integrand and code. – MarcoB Dec 30 '19 at 20:28
• @MarcoB In principle, it is a valid question if we take it as written: the OP simply asks for an explanation of the actual meaning of the NIntegrate::eincr message. – Alexey Popkov Dec 31 '19 at 9:14
• @Alexey The meaning of the error is straightforward and it is explained in the error message. The actual question is "Should I believe in NIntegrate more than Cuba in this case? Or both are just bad results?"; unfortunately that question is not easily answerable without more details. – MarcoB Dec 31 '19 at 18:41
• Following is an answer I began before the question was closed. I'll post it as an answer here (and proofread it!), if the question is reopened. github.com/mroge02/math315/blob/master/MaxErrorIncreases.md -- arguably not worth it, but it's a good tutorial, ihmo. It's a response to "Just want to know the meaning of 'global error of the strategy GlobalAdaptive has increased more than X times'." (I explain only in terms of 1D integrals, but it's easier to understand and illustrate.) – Michael E2 Jan 1 at 6:14

## Introduction

I'll give a generic answer to a generic question (and mostly ignore Cuba).

What is supposed to happen in the global-adaptive strategy is that each recursive refinement, which subdivides the interval with the largest integration error, should reduce the error estimate. But it doesn't always happen. Why not? Let's consider a few common, single-variable cases:

• The integrand is analytic or at least sufficiently smooth (the derivatives exist and are continuous to a higher order that the integration rule).
• The integrand has a singularity in the interval of integration.
• The integral is zero or near zero.
• The integrand is highly oscillatory.

Before we get started, let me mention that the integral and error estimate are based on the evaluation of the integrand $$f(x)$$ at a discrete set of nodes $$x_1,\dots,x_n$$. Note that $$f(x) + (x-x_1)\cdots(x-x_n)\,g(x)$$ has the same function values at the nodes as $$f$$ and therefore has the same integral and error estimate as $$f$$. Since $$g(x)$$ is arbitrary, the actual error (not the computed error estimate) can be arbitrary. In practice it is rare that you get this sort of aliasing problem except for a few of the subdivisions, since the nodes change when a subinterval is divided. (The bigger danger here is that the error is falsely estimated as close to zero and the subinterval is never divided. But that has nothing to do with the NIntegrate::eincr warning.)

## Four cases

### Smooth integrand

In interpolatory function approximation, which is used in most integration rules, there is a pre-convergent phase followed by a convergent phase.

Example. A three-point interpolatory rule (using a quadratic polynomial interpolant) cannot make a good approximation to a function that has five local extrema (or inflection points). Either you have to use a higher-order rule (use more points) or subdivide the interval so that in each subinterval the function has at most one extremum; and if an inflection point cannot be avoided, make the interval very small so that the bad approximation will have little effect. (Those of you who know Simpson's rule know that because of its symmetry when integrated, we get a superconvergent rule with degree of precision 3 that can integrate cubic inflection points exactly. But an asymmetric 3-point rule won't have this superconvergence, so bear with me.) The main point is that the global-adaptive method has to break up the interval into small enough pieces that the integration rule starts to make good approximations of the integrand or its integral. After that point, the method will enter the convergent phase.

During the preconvergent phase, one can expect the error estimate to increase sometimes. For a highly oscillatory integrand, it might increase often. If it increases too many times (for NIntegrate), then raising "MaxErrorIncreases" would theoretically fix the problem (and practically would, too, if the number of increases were not so great as to take too much time). Another approach is to increase the order of the interpolatory rule using the "Points" suboption (see the docs; beware that Newton-Cotes becomes less stable as the order increases); raising the order can shorten the preconvergent phase.

### Singularity

A singularity often results in a NIntegrate::slwcon warning, and sometimes in a NIntegrate::eincr one, too. Singularities tend to yield larger error estimates. As an interval is subdivided, the effect of the singularity might be come more noticeable. Another way of saying it is that singularities can extend the preconvergent phase and error increases can accumulate. The symbolic preprocessing stage of NIntegrate often identifies the singularities often locates the singularities, but they can be user-specified with the "Exclusions" option or by listing them in the integration-interval iterator, e.g., {x, a, s1, s2,..., sk, b}.

### Zero integrals

Convergence to the desired PrecisionGoal can be particularly troublesome for integrals that are equal to zero or very close to zero (compared to the amplitude of the function, the working precision, and the PrecisionGoal). Recall that a precision goal $$p$$ is met by the error being smaller than the result by a factor of $$10^{-p}$$. Potentially in an near-zero integral, round-off error will dominate the truncation error of the integration method, and the error will jitter about causing many error increases. This last problem is handled by setting AccuracyGoal to a finite value $$a$$ such that you are comfortable that an error estimate less $$10^{-a}$$ should be accepted as close enough to zero to be considered zero and a further increase of precision should not be sought.

### Oscillatory integrand

This is actually a special subcase of a smooth integrand. Oscillatory integrals usually require a special method, since they have infinitely many extrema, and you can't just break up the interval into pieces in which the number of extrema is somewhat less than the order of the integration rule. The preconvergent phase could be arbitrarily long, depending on how fast the amplitude decays (which could depend on a transformation), during which error increases may be numerous.

## Examples

I've picked some contrived, simple-ish examples that can be investigated analytically. In some cases, I force a bad method to illustrate how eincr works.

### Smooth integrand, long preconvergent phase

nn = 2406;
ff[x_?NumericQ] := ChebyshevT[nn, x]^2;
NIntegrate[ff[x], {x, -1, 1}, MaxRecursion -> 16,
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]
% - Integrate[Cos[nn t]^2 Sin[t], {t, 0, Pi}] (* exact error *)


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: ...highly oscillatory integrand....

NIntegrate::eincr: ...NIntegrate obtained 0.9294557079749661 and 0.0013688527738452687 for the integral and error estimates.

(*
0.929456
-0.0705442
*)


Note the exact error is quite a bit larger than the error estimate. This can be fixed with "MaxErrorIncreases":

NIntegrate[ff[x], {x, -1, 1}, MaxRecursion -> 16,
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule",
"MaxErrorIncreases" -> 1000}]
% - Integrate[Cos[nn t]^2 Sin[t], {t, 0, Pi}]  (* exact error *)


NIntegrate::slwcon: ....

(*
0.9999999565975315
3.32068*10^-13
*)


The problem can also be handled by pre-subdividing the interval. The error estimate is not computed for this, and so it does not count in the number of error increases. It also does not count in the heuristics that lead to the NIntegrate:slwcon warning, and sometimes it prevents that message. One way is to use MinRecursion:

NIntegrate[ff[x], {x, -1, 1}, MinRecursion -> 10, MaxRecursion -> 16,
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]
% - Integrate[Cos[nn t]^2 Sin[t], {t, 0, Pi}]
(*
1.
1.62093*10^-14
*)


Another way is to subdivide the interval manually at, say, the zeros of ff[x]:

NIntegrate[ff[x],
Evaluate@Flatten@{x, Join[{-1}, -Cos[Pi ( Range[nn] - 1/2)/nn], {1}]},
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]
% - Integrate[Cos[nn t]^2 Sin[t], {t, 0, Pi}]
(*
0.9999999568134225
-6.66134*10^-16
*)


### Singularities

The following integrand has weak singularities where Sin[x] changes sign (the derivative of the integrand is undefined).

integrand[x_?NumericQ] := Abs[Sin[x]]; (* prevents symbolic preprocessing *)
NIntegrate[integrand[x], {x, 0, 800 Pi}, MaxRecursion -> 20]

> NIntegrate::eincr: ...NIntegrate obtained 1600.083569906413 and 0.18856172376838135 for the integral and error estimates.

(*  1600.08  *)


This can be handled by pre-subdividing the interval. The error estimate is not computed for this, and so it does not count in the number of error increases. (Note this approach can be used when there is a long preconvergent phase, too.)

NIntegrate[integrand[x], {x, 0, 800 Pi},
MinRecursion -> 10, MaxRecursion -> 20]
(*  1600.  *)

NIntegrate[integrand[x], Evaluate@Flatten@{x, Range[0, 800] Pi}]
(*  1600.  *)


### Integral near zero

NIntegrate[Sin[x], {x, 0, 2 Pi}, MaxRecursion -> 16]


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: ...value of the integration is 0....

NIntegrate::eincr: ...NIntegrate obtained 4.631711680858075*^-16 and 5.169593219712382*^-16 for the integral and error estimates.

(*  4.63171*10^-16  *)


A fix is to set an AccuracyGoal a bit larger than the round-off error:

NIntegrate[Sin[x], {x, 0, 2 Pi}, MaxRecursion -> 16,
AccuracyGoal -> 15]
(*  -2.6507*10^-16  *)


### Smooth oscillatory integrand, nonzero integral

This example has an unexpected "coincidence," which took me off guard. But it shows that when numerical methods go bad, you should be open and separate what you actually know from what usually happens. The interval comprises many, many periods of the integrand. Predictably, during the first several subdivisions, the error bounces around. Since the interval is so long, it takes many subdivisions during which too many error increases occur.

NIntegrate[Sin[x]^2, {x, 0, 2^12 Pi}, MaxRecursion -> 11,
Method -> "GaussKronrodRule"]
% - 2^12 Pi/2. (* actual error *)


NIntegrate::eincr: ...NIntegrate obtained 7328.140790022457and 295.19159332164276 for the integral and error estimates.

(*
7328.14
894.159
*)


A slight change in the interval and we get a very good estimate of the integral automatically. This turns out to be a surprising accident. It computes the correct result with just one subdivision, which you can check by adding MaxRecursion -> 1.

NIntegrate[Sin[x]^2, {x, 0, (2^12 - 1) Pi},
Method -> "GaussKronrodRule"]
% - (2^12 - 1) Pi/2. (* actual error *)
(*
6432.41
-5.63887*10^-11
*)


Change it slightly again, just by adding some preliminary subdivision, and we're in trouble again.

NIntegrate[Sin[x]^2, {x, 0, (2^12 - 1) Pi}, MinRecursion -> 2,
Method -> "GaussKronrodRule"]
% - (2^12 - 1) Pi/2. (* actual error *)


NIntegrate::eincr: ...NIntegrate obtained 6377.785070697375 and 143.07090946442491 for the integral and error estimates.

(*
6377.79
-54.6259
*)


Just to illustrate that the Gauss-Kronrod rule I chose above is contrived to produce the NIntegrate::eincr problem: The Levin rule is the automatic rule chosen by NIntegrate, and it produces a good result, with much more than the eight digits of precision sought with the default PrecisionGoal.

NIntegrate[Sin[x]^2, {x, 0, 2^12 Pi}, Method -> "LevinRule"]
% - 2^12 Pi/2. (* actual error *)
*&
6433.98
6.89033*10^-9
*)


No doubt you're curious what is going on. Here is the coincidence: It happens that a (pretty useless) property of the Gauss-Kronrod rule is that it calculates the integral of $$\sin^2 x$$ over an interval of the form $$[m \pi/2, n \pi/2]$$, where $$m$$ and $$n$$ are positive integers such that $$m+n$$ is odd, exactly and estimates the error to be zero (we don't get zero exactly because of round-off error; but run it with WorkingPrecision -> 16 and you do get zero). In the second integral, the first subdivision creates two subintervals of that form, $$[0, (2^{12}-1)\pi/2]$$ and $$[(2^{12}-1)\pi/2, 2^{13} \pi/2]$$, and the integral is computed exactly over each of them. In the third integral, the interval is subdivided twice before any integration is computed. This bypasses the special form intervals, and the Gauss-Kronrod rule is no longer exact. NIntegrate proceeds as in the first integral, and we get an eincr message.

## Visualization

NIntegrate comes with an undocumented tool, the option IntegrationMonitor. With it, we can see the error increases (and could even show which subdivisions cause each one, but I'll omit that).

The following is a function that uses IntegrationMonitor to plot the error after each subdivision. You can use it to illustrate any of the examples above. You can also use it on any NIntegrate command whose method or rule utilizes IntegrationMonitor (not all do).

ClearAll[errorPlot];
SetAttributes[errorPlot, HoldAll];
errorPlot[nint_NIntegrate, plotopts___?OptionQ] := Block[{integral, errors},
{integral, errors} = Reap[
Hold[nint] /.
Hold[Verbatim[NIntegrate][args___]] :> NIntegrate[args,
IntegrationMonitor :> (Sow[Total@Through[#@"Error"]] &)]
];
errors = Flatten@errors/integral;
ListLinePlot[errors // RealExponent, PlotRange -> All,
Frame -> True, FrameLabel -> {"subdivisions", "log error"},
PlotLabel ->
Row[{UnitStep@Differences@errors // Total, " error increases"}],
DataRange -> {0, Length@errors - 1},
plotopts]
]


We'll apply it to the last example. First a shorter interval, so we can see and count the error increases easily. After around 16 subdivisions, the preconvergent phase ends. One can check that in this (very symmetric) integral, we have 16 equal subintervals. To see the convergent phase, each of these has to be subdivided. After 32 subdivisions the error decreases dramatically. It will do it again after 64 subdivisions. (There is an obvious blip at subdivision 31, but I did not investigate it. You can use IntegrationMonitor to do so, if curious.)

errorPlot[NIntegrate[Sin[x]^2, {x, 0, 2^5 Pi},
Method -> "GaussKronrodRule"]]


In the main example, we see that we're still in the preconvergent phase after 800+ subdivisions.

errorPlot[NIntegrate[Sin[x]^2, {x, 0, 2^12 Pi},
Method -> "GaussKronrodRule", MaxRecursion -> 11]]


Going further, we see that the preconvergent phase ends after about 2000 subdivisions (or 2^11), and a big leap in convergence happens after another 2000 steps. A second leap occurs after another 4000 steps. (Please keep in mind this doubling of the number of subdivisions comes from the symmetry of the integral and is not at all typical.)

errorPlot[NIntegrate[Sin[x]^2, {x, 0, 2^12 Pi},
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule",
"MaxErrorIncreases" -> 1000}, MaxRecursion -> 16]]
`

• This is a nice article. – Artes Mar 6 at 22:34