Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I need to numerically integrate a highly oscillatory function over the semi-infinite domain $(0,\infty)$:

$$\int_0^\infty \frac{\sin^2(x) \sin^2(1000 x)}{x^{5/2}}\mathrm dx$$

Since the Levin rule (which was recently added to Mathematica, starting in version 8) was developed specifically for oscillatory integrals such as this, I thought I'd try it:

ans = NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity}, 
                 Method -> {"LevinRule"}, PrecisionGoal -> 8, MaxRecursion -> 30]

Using an exact solution for this integral, I can confirm the relative accuracy of the Mathematica result is $1 \times 10^{-11}$, and moreover the calculation is very quick. At first, this led me to believe that Levin's method works great for this problem, but...

It turns out that Mathematica must be automatically switching to non-oscillatory rule behind the scenes, because forcing it not to do so gives a very poor result:

ans = NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity},  
                 Method -> {"LevinRule", "MethodSwitching" -> False},
                 PrecisionGoal -> 8, MaxRecursion -> 30]

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 30 recursive bisections in x near {x} = {0.}. NIntegrate obtained -3497.5 and 3510.0321785369356` for the integral and error estimates. >>

Is there any way to find out which alternative non-oscillatory rule Mathematica is automatically selecting? I've tried to guess which rule is being used by manually specifying few rules but the results I've obtained with other rules are inaccurate, slow, or both:

ans = NIntegrate[(Sin[x])^2 (Sin[1000 x])^2/x^(5/2), {x, 0, Infinity}, 
                 Method -> "ClenshawCurtisRule", AccuracyGoal -> 8, MaxRecursion -> 30]

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 400 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.000013202052151832003` and 1.0480362255168103`*^-6 for the integral and error estimates. >>

I'd like to know what rule Mathematica is using so that I can try adjusting options for the best performance possible. I need to calculate this integral several hundred thousand times, as the innermost integral of a nested double integration. Furthermore, when it comes to publishing my results, I would like to be able to state the integration strategy that was actually used, rather than "Mathematica knew how to handle it".

share|improve this question
Why not use, say, "ExtrapolatingOscillatory" or "DoubleExponentialOscillatory" instead, since your oscillatory integral is over $[0,\infty)$ anyway? – J. M. Jun 4 '13 at 17:20
Thanks for the quick response! There are a few reasons I didn't use "ExtrapolatingOscillatory" or "DoubleExponentialOscillatory" from the start. First, I will need to integrate this function over both finite intervals [0,a] as well as semi-infinite intervals [a,infinity], with a>0. Second, I had previously implemented both of these rules from scratch in IDL and found that both require too much computation time when the frequency (1000 in the example I show here) increases. For example, to compute – user7885 Jun 4 '13 at 17:40
BTW: in fact, the default method used by NIntegrate[] for "nice" integrals is "GlobalAdaptive"; as you've noticed, it is the method being used instead of "LevinRule", somewhat sneakily... – J. M. Jun 4 '13 at 17:46
Sorry--my previous comment was truncated. I didn't opt for the "ExtrapolatingOscillatory" or "DoubleExponentialOscillatory" rules because I will need to integrate this function over both finite intervals [0,a] as well as semi-infinite intervals [a,infinity]. Also, computation time for these rules should scale with frequency (w=1000 in the example I show) but not so for Levin. Nevertheless, having just tried it, I see that both the "ExtrapolatingOscillatory" and "DoubleExponentialOscillatory" rules in Mathematica do give fast and accurate results even for w=100,000. – user7885 Jun 4 '13 at 17:52
In particular, you might be interested in the ExcludedForms option, if you choose the strategy of only considering the more oscillatory Sin[1000 x], while ignoring Sin[x]. If you take a look at the advanced documentation, there is also a pointer to the undocumented function NIntegrate`LevinIntegrandReduce[], which might help you in exploring how to properly set options. – J. M. Jun 4 '13 at 18:49

2 Answers 2

As you may know, the backup non-oscillatory rule is controlled by the Method sub-option of `"LevinRule", documented here. You can use it like this:

In[60]:= NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity}, 
 Method -> {"LevinRule", 
   Method -> {"GaussKronrodRule", "Points" -> 11}}]

Out[60]= 1.16762

With this, you can try to tune for performance.

In your example, the default setting Method -> Automatic is effectively equivalent to Method -> {"GaussKronrodRule", "Points" -> 5}. There is no documented method to determine this, but you can more or less verify it by doing something naughty like:

In[61]:= Block[{NIntegrate`GaussKronrodRuleData = (Print[{##}]; 
     Abort[]) &}, 
 NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity}, 
  Method -> {"LevinRule"}, PrecisionGoal -> 8, MaxRecursion -> 30]]

During evaluation of In[61]:= {5,MachinePrecision}

Out[61]= $Aborted
share|improve this answer

This question comes up often enough. See this discussion at : Integration method used in NIntegrate , and the notebook Finding the applied NIntegrate methods attached to my second response in the discussion.

That notebook contains examples of usage of the undocumented function NIntegrateSamplingPoints and NIntegrate's option IntegrationMonitor.

The integral in the question

For the integral in the question with NIntegrateSamplingPoints we get kind of a boring picture because of the infinite region. Taking logs of the sampling points might be more informative:


gr = NIntegrateSamplingPoints@
   NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity}, 
    Method -> {"LevinRule", 
      Method -> {"GaussKronrodRule", "Points" -> 11}}];

Graphics[gr[[1]] /. 
  Point[{x_?NumericQ, y_?NumericQ}] :> 
   Point[{Log[10, x + 10^-12], y}], Frame -> True, 
 FrameLabel -> {"lg(sampling points)", "evaluation order"}, AspectRatio -> 1/1.5]

enter image description here

The plot shows the order of evaluation of the sampling points.

Using IntegrationMonitor we can see the application of the integrand over the regions derived with the LevinRule method:

t = Reap[NIntegrate[Sin[x]^2 Sin[1000 x]^2/x^(5/2), {x, 0, Infinity}, 
    Method -> {"LevinRule", 
      Method -> {"GaussKronrodRule", "Points" -> 11}}, 
    PrecisionGoal -> 1.5, 
    IntegrationMonitor -> (Sow[
        Map[{#1@"Integrand", #1@"Boundaries", #1@"Integral", #1@
            "Error"} &, #1]] &)]];
res = t[[1]];
t = t[[2, 1]];

enter image description here


(I have had the code below for some time, but I have been hesitant to share it because of several concerns. It is somewhat hard to interpret its results and coming up with it does require some internal knowledge of NIntegrate's development. After dedicated discussions with/about the NIntegrate method tracing code at WTC 2015 it seems that it is better to show and describe it.)

We can trace NIntegrate's method initialization by manipulating the top level initialization function implementations. The basic idea is to take the down values and up values of the initialization functions of NIntegrate's methods that have the form

Block[{v___}, b_CompoundExpression]

and replace them with

Block[{res = Block[{v}, b]}, Print[res]; res]

When an NIntegrate command is executed we will see a printed trace of the methods that are initialized.

Here is the tracing code:

symbNames = Names["NIntegrate`*"];
symbNames = 
     symbNames, (__ ~~ "Rule") | (__ ~~ 
        "Global" | "Local" | "MonteCarlo" | "Principal" | "Levin" | 
         "Osc" ~~ ___)]], "NIntegrate`AutomaticStrategy"];
symbs = ToExpression[#] & /@ symbNames;
dvs = DownValues /@ symbs;
uvs = UpValues /@ symbs;
Unprotect /@ symbs;
dvsNew = MapThread[
   With[{s = #2},
     DownValues[s] = 
         a_ :> b___] :> (a :> (Print["DownValue call for: ", 
            Style[s, Red]]; b))]] &, {dvs, symbs, symbNames}];
uvsNew = MapThread[
   With[{s = #2},
     UpValues[s] =
       HoldPattern[Block[vars_, CompoundExpression[b___]]] :>
        Block[{res = Block[vars, CompoundExpression[b]]}, 
         Print["UpValue call for: ", Style[s, Blue], 
          Style[" ::\n", Blue], res]; res]
     ] &, {uvs, symbs, symbNames}];

In that code using Pick I have reduced the number of NIntegrate context symbols being traced. Of course, if it is desired, the down/up values of the full list of the NIntegrate context symbols can be manipulated for tracing.

Let us look at tracing examples.

Here is a numerical integration computation with an automatically selected method:

enter image description here

The basic objects of NIntegrate are integration regions. Each region has its own integration function and integration rule. In order to interpret the printed trace it is helpful to know that NIntegate uses the software design patterns Strategy, Composite, Decorator, and others.

NIntegrate's symbolic pre-processors are using Decorator, and we can see in the trace that the outcome of AutomaticStrategy are pre-processor symbols wrapped around the method GlobalAdaptive. The method global adaptive uses a Gauss-Kronrod rule, which after its initialization is treated as a general rule. (I.e. a list of abscissas, integral estimate weights, and approximation error weights.)

The method GlobalAdaptive is going to be directly used if symbolic processing is prevented:

enter image description here

Here is a numerical integration computation with a specifically selected pre-processing method:

enter image description here

share|improve this answer
It seems you know how to work with Experimental`CreateNumericalFunction. Please consider answering the dedicated question or improving my CW answer in that thread. It is currently unclear how to use the fourth argument of CreateNumericalFunction as well as options "ErrorReturn", "SampleArgument" and Message. As you can see, the question is highly appreciated by the community. – Alexey Popkov Oct 10 at 12:14
@AlexeyPopkov I am sorry I don't think I can contribute more to that post. (And my knowledge on that function seems to be fairly outdated...) – Anton Antonov Oct 11 at 0:28
It's too bad I can't give another +1. Thanks for sharing! – J. M. Nov 4 at 4:39
@J.M. Thanks, that is nice to hear! – Anton Antonov Nov 4 at 12:14

Your Answer


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.