5
$\begingroup$

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?

$\endgroup$
1
3
$\begingroup$

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:

with Sin

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

with bessel

$\endgroup$
5
  • 1
    $\begingroup$ 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...) $\endgroup$ – J. M.'s ennui Jun 5 '20 at 17:23
  • $\begingroup$ 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. $\endgroup$ – Jolyon Jun 5 '20 at 17:43
  • 1
    $\begingroup$ @Jolyon, you might want to look at how to use "LevinRule" if you expect to be doing a lot of integrals with spherical Bessels. $\endgroup$ – J. M.'s ennui Jun 5 '20 at 17:53
  • $\begingroup$ Use FunctionExpand to simplify the integrand before integration. NIntegrate[ 100 SphericalBesselJ[0, 100 k]/k^4 // FunctionExpand, {k, 100, 250}] $\endgroup$ – Bob Hanlon Jun 5 '20 at 19:23
  • $\begingroup$ @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. $\endgroup$ – J. M.'s ennui Jun 6 '20 at 0:17

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.