# Determining which rule NIntegrate selects automatically

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".

-
Why not use, say, "ExtrapolatingOscillatory" or "DoubleExponentialOscillatory" instead, since your oscillatory integral is over $[0,\infty)$ anyway? –  Ｊ. Ｍ. 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... –  Ｊ. Ｍ. 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 NIntegrateLevinIntegrandReduce[], which might help you in exploring how to properly set options. –  Ｊ. Ｍ. Jun 4 '13 at 18:49

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[{NIntegrateGaussKronrodRuleData = (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

-