11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Look e.g. at this answer How to numerically integrate this integral where you can find an automatic error estimate and two recipes how to deal with this problematic issue. $\endgroup$ – Artes Feb 19 '15 at 18:16
  • $\begingroup$ I'd hope that there is some 'real' way to get the estimate, but you could force a message containing the error estimate by reducing the value of MaxPoints and then catching the value of the message NIntegrate::maxp. $\endgroup$ – 2012rcampion Feb 19 '15 at 18:27
  • $\begingroup$ See tutorial/NIntegrateIntegrationRules (search for error weights). See also the rule data functions ?NIntegrate*RuleData`, which return the error weights, among other things. I used such functions in my answer to 61238 to calculate the error estimate. $\endgroup$ – Michael E2 Feb 19 '15 at 21:37
  • $\begingroup$ @Artes: As far as I can tell that will not allow me to get me the error estimate when no error messages are generated. $\endgroup$ – mmeent Feb 20 '15 at 9:15
  • $\begingroup$ @2012rcampion : In principle that is a good idea. However, that will generally cause NIntegrate to stop before having reached the desired accuracy. If one were doing a single numerical integral, one could fiddle with the settings to generate the error in this way. However, in the current application I am think of, there will be tens of thousands of NIntegrates performed, and I want to generate error bars on their results. $\endgroup$ – mmeent Feb 20 '15 at 9:20
11
$\begingroup$

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
*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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