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?

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
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

Improve this answer
$\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 eventual burnout ♦
    Jun 5, 2020 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, 2020 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 eventual burnout ♦
    Jun 5, 2020 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, 2020 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 eventual burnout ♦
    Jun 6, 2020 at 0:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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