Obtaining an NIntegrate error estimate

Question

Is there a way to extract the error that Mathematica estimates when calculating a numerical integral using NIntegrate?

Internally Mathematica must keep track of this error, because it is used to determine if the PrecisionGoal has been met.

The reason I want to extract this information is that some integration strategies (e.g. "Trapezoidal") can significantly overshoot the PrecisionGoal target. In cases where this happens it would be very useful to include this in any estimated error bars on the result.

edit

For clarification lets consider example.

I have a complicated integrand f[x] (that is expensive to compute). I know a couple of things about this integrand:

1) It is 2π periodic. 2) It is C-infinite smooth.

These two facts mean that f[x] has a Fourier series whose coefficients decay exponentially.

This in turn means that a trapezoidal integration strategy also converges exponentially. Hence I integrate with:

NIntegrate[
 f[x],
 {x,0,2Pi},
 Method-> {"Trapezoidal", "SymbolicProcessing"->0},
 PrecisionGoal -> n
]

Where n is my desired precision. This works well.

Now I want to estimate the error bound of the integration. More I want to do this without any further evaluations of f[x].

Because (by assumption) of exponential convergence of the integrand, the "Trapezoidal" strategy will double the precision of the result at each integration step (once it is in the tail). Consequently, the (estimated) precision of the final answer is somewhere between n and 2n. By guessing an error bound equal to the precision goal we my be massively overestimating the actual error. (Which is not good if the result is to be used in a later data analysis step.)

This integral is part of a loop of a much longer code. A typical run contains up to 10^5 of these integrals. Hence, fiddling with settings to coax an error report for a single integral is not really an option. Somewhere in its internals Mathematica is calculating this error estimate, hence it must be possible to extract it. If only we knew the name of the internal variable that is used for the error estimate for the trapezoidal strategy.

Answer 1

From Anton Antonov's answer to Determining which rule NIntegrate selects automatically, we learn we can get the error estimate(s) using the undocumented option IntegrationMonitor.

Here's test example adapted from the tutorial on the "Trapezoidal" strategy.

f[x_] := 1/π Cos[80 Sin[x] - x];

exact = Integrate[f[x], {x, 0, 2 Pi}]
(*  2 BesselJ[1, 80]  *)

res = NIntegrate[f[x], {x, 0, 2 Pi}, 
  Method -> {"Trapezoidal", "SymbolicProcessing" -> 0},
  IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)]
Total@errors
Abs[res - exact]
(* 
  -0.112115          <-- integral value
  2.67841*10^-15     <-- error estimate by NIntegrate
  4.996*10^-16       <-- actual error
*)

Note the errors are given in list of error estimates, one for each subregion of integration. The "Trapezoidal" generates only one region. But in general one needs to total these errors. For instance, here the default "GlobalAdaptive" strategy generates 175 subregions:

NIntegrate[f[x], {x, 0, 2 Pi},
 PrecisionGoal -> 8,
 IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)]
Length@errors
Total@errors
(*
  -0.112115          <-- integral value
  175                <-- length of errors = number of subregions
  1.04099*10^-9      <-- overall error estimate
*)