Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am trying to evaluate a highly oscillatory integral using NIntegrate. I fear that due to limited resources (time and/or memory), I will not be able to evaluate the integral to the desired precision. Thus, I would like to programmatically access the error estimates that are e.g. reported by the messages NIntegrate::maxp, NIntegrate::ncvb, or NIntegrate::eincr. I could not find an option of NIntegrate that would directly make these error estimates available. However, given that I have to evaluate a multitude of integrals, it is impractical to obtain the errors from the warnings by hand.

The following example generates the NIntegrate::maxp message (obviously this very integral has an analytical solution):

NIntegrate[Sin[x]/Sqrt[x], {x, 0, 100}, Method -> "MonteCarlo",PrecisionGoal -> 6]

NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 1.1787733508261242and 0.07678430788995934 for the integral and error estimates.

How to get (if necessary, extract) the error estimate (0.07678430788995934`)?

Remark: The example from the help of NIntegrate::eincr, i.e. ref/message/NIntegrate/eincr, does not produce the expected message in version 8.0; unfortunate my integrals still do.

share|improve this question
2  
"...evaluate a highly oscillatory integral using NIntegrate[]." - then, why "MonteCarlo"? There's "DoubleExponential" or "ClenshawCurtisOscillatoryRule" which you could have used... unless your actual integrals are in fact multidimensional, and you've just grossly oversimplified. – J. M. Nov 13 '12 at 11:39
2  
Yes. I am grossly simplifying and the actual integral is multidimensional (4D). In particular, I have just chosen the method since it generates one of the messages in question. It turned out that generating the NIntegrate::eincr message with a simple 1D integral is unexpectedly (given my mathematical naivety) difficult, i.e. NIntegrate is very robust (see my remark). For the actual integral, Method->”MonteCarlo” is in fact my best bet. – dan Nov 13 '12 at 12:54
up vote 12 down vote accepted

It seems to me that there's a better approach, but one way is to define your own DownValue for this particular message. For example:

Unprotect[Message];
Message[NIntegrate::maxp, its_, int_, err_] := Sow[err]

Then

NIntegrate[Sin[x]/Sqrt[x], {x, 0, 100}, 
  Method -> "MonteCarlo", PrecisionGoal -> 6] // Reap

(* Out: {1.07721, {{0.0761274}}} *)
share|improve this answer
1  
Excellent. For NIntegrate::eincr the same approach applies. For NIntegrate::ncvb, Message[NIntegrate::ncvb, nr_, var_, varlist_, at_, int_, err_] := Sow[err] seems to do the job. In general, Message[NIntegrate::ncvb | NIntegrate::eincr | NIntegrate::maxp, l___] := Sow[Last@{l}] could probably be used. – dan Nov 13 '12 at 13:25
1  
In order to retain the original message, one could probably use something along the lines of: Block[{original = False}, Unprotect[Message]; Message[NIntegrate::maxp, l___] /; Not[original] := (Sow[Last@{l}]; original = True; Message[NIntegrate::maxp, l]); NIntegrate[Sin[x]/Sqrt[x], {x, 0, 100}, Method -> "MonteCarlo", PrecisionGoal -> 6] // Reap ] – dan Nov 13 '12 at 13:33

From the description of the question it seems to me that using the (undocumented) option IntegrationMonitor to obtain integration intervals and estimates might be very useful.

Here is an example:

t = Reap[NIntegrate[Sin[x]/Sqrt[x], {x, 0, 100}, PrecisionGoal -> 6, 
    Method -> "MonteCarlo", 
    IntegrationMonitor -> (Sow[
        Map[{#1@"Boundaries", #1@"Integral", #1@"Error"} &, #1]] &)]];
res = t[[1]];
t = t[[2, 1]];
Take[t, -4]

enter image description here

More examples and explanations about the use of IntegrationMonitor can be found in the notebook "Finding the applied NIntegrate methods.nb" attached to the community.wolfram.com discussion "Integration method used in NIntegrate".

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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