# Sine vs Sinc vs SphericalBesselJ in NIntegrate

I'm evaluating an oscillatory integral numerically, and ran into a weirdness with NIntegrate, which I've boiled down to a simple case for this question.

Consider NIntegrate[Sin[100 k]/k^5, {k, 100, 250}], which evaluates to -9.59482*10^-13. Mathematica can handle this integral analytically, writing the solution in terms of Sine integrals. Doing numerical evaluation of the result yields the same answer to 13 digits (N[Integrate[Sin[100 k]/k^5, {k, 100, 250}], 15]).

Now, consider writing the Sin function instead as a Sinc function: NIntegrate[100 Sinc[100 k]/k^4 , {k, 100, 250}]. Same result.

Finally, consider writing the Sinc function as a SphericalBesselJ function: NIntegrate[100 SphericalBesselJ[0, 100 k]/k^4 , {k, 100, 250}]. Note that SphericalBesselJ[0,x] == Sinc[x]. Here, Mathematica complains bitterly about doing this integral numerically, before offering the answer of 6.43722*10^-11.

Why isn't Mathematica choosing the same method to numerically integrate SphericalBesselJ as it is for Sin and Sinc?

NIntegrate may be hardwired to "recognizes" structure in Sin and Sinc and automatically select better methods, but its heuristics are not smart enough to do that with SphericalBesselJ. I would not have known either :-).

Perhaps in further support of that the following executes with no trouble and returns the same results as the others:

NIntegrate[100 SphericalBesselJ[0, 100 k]/k^4, {k, 100, 250}, Method -> "LevinRule"]
(* Out: 9.59482*10^-13 *)


Edit

To go qualitatively a little further, we can use the beautiful spelunking code provided by Anton Antonov in his answer to Determining which rule NIntegrate selects automatically. Once you execute his tracing code from that answer, then the following calls return a lot more info:

Note that the above includes a call to NIntegrateClenshawCurtisOscillatoryRule. This rule is NOT invoked in the SphericalBesselJ call:

• As an additional note: "ClenshawCurtisOscillatoryRule" (see e.g. this paper) is specifically designed for finite oscillatory integrals with explicit sine or cosine factors, so of course NIntegrate[] would not consider using it for an integrand with SphericalBesselJ[], even if it is equivalent. (Perhaps Anton might have more to say about this matter...) Jun 5, 2020 at 17:23
• Well, isn't that a curiosity. Guess I should be replacing SphericalBesselJ with Sinc explicitly wherever I can to help out with the method selection, or suggesting an integration method for Mathematica. Jun 5, 2020 at 17:43
• @Jolyon, you might want to look at how to use "LevinRule" if you expect to be doing a lot of integrals with spherical Bessels. Jun 5, 2020 at 17:53
• Use FunctionExpand to simplify the integrand before integration. NIntegrate[ 100 SphericalBesselJ[0, 100 k]/k^4 // FunctionExpand, {k, 100, 250}]` Jun 5, 2020 at 19:23
• @Bob, here it works OK, but in general, I would not recommend converting spherical Bessels to their trigonometric form. See e.g. this and this. Jun 6, 2020 at 0:17