Why NIntegrate is badly-behaved on $J_{\frac{9}{2}}(x)$ by default?

Question

A friend of mine showed me this example:

Plot[BesselJ[9/2, x], {x, 0, 1}, 
 PlotLabel -> Style["The integrand seems to be simple", 14]]

Integrate[BesselJ[9/2, x], {x, 0, 1}] // N
(* 0.000148473 <- This is the correct result *)

NIntegrate[BesselJ[9/2, x], {x, 0, 1}]
(*-1.26625*10^170 *)

It's not hard to fix the code:

(* Solution 1 *)
NIntegrate[BesselJ[9/2, x], {x, 0, 1}, WorkingPrecision -> 16]
(* 0.0001484729674125616 *)

(* Solution 2 *)
Clear@f
f[x_?NumericQ] := BesselJ[9/2, x]
NIntegrate[f@x, {x, 0, 1}]
(* 0.000148473 *)

But I wonder why NIntegrate fails in such a wild way by default?

Seeing the Solution 2 above, it's natural to guess that NIntegrate has done some improper symbolic processing, but this seems not to be true:

NIntegrate[BesselJ[9/2, x], {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}]
(* -1.26625*10^170 *)

Answer 1

Because NIntegrate evaluates the integrands before starting the actual integration, in some cases (like this one) it is better to define the integrand function F with the signature F[x_?NumericQ].

BF[n_?NumericQ, x_?NumericQ] := BesselJ[n, x]

NIntegrate[BF[9/2, x], {x, 0, 1}]

(* 0.000148473 *)

Integrate[BesselJ[9/2, x], {x, 0, 1}]
%% // N

(* Sqrt[2/\[Pi]] (30 Cos[1] + 
   Sqrt[2 \[Pi]] FresnelS[Sqrt[2/\[Pi]]] - 20 Sin[1]) *)

(* 0.000148473 *)

When expanded

BesselJ[9/2, x] // Expand

(* -((105 Sqrt[2/\[Pi]] Cos[x])/x^(7/2)) + (
 10 Sqrt[2/\[Pi]] Cos[x])/x^(3/2) + (105 Sqrt[2/\[Pi]] Sin[x])/x^(
 9/2) - (45 Sqrt[2/\[Pi]] Sin[x])/x^(5/2) + (
 Sqrt[2/\[Pi]] Sin[x])/Sqrt[x] *)

the oscillatory nature of the integrand terms near zero confuses NIntegrate's algorithms.

(I have answered very similar questions in different forums over the years.)

A similar discussion/answer is J.M.'s answer of "NIntegrate and Integrate of a Spherical Bessel function".

And as J.M. noted in a comment :

"SymbolicProcessing" -> 0 does absolutely nothing because BesselJ[] is already auto-expanded in the half-integer case, and the resulting expression is numerically iffy.

(We mention this because of the attempt to remedy the situation with "SymbolicProcessing" in the question.)

Answer 2

I tried to make it a comment and then can't control the words. I think the problem is arising from the fact that the value of the function and its derivative is too small to near x=0.

Plot[BesselJ[9/2, x], {x, 0, .01}] 
Plot[Evaluate[D[BesselJ[9/2, x], x]], {x, 0, 0.01}]

As you can see the derivative is highly oscillatory at this small range.

When you are using Integrate it is doing a symbolic evaluation which keeps this violent things out of the picture. When you are going with NIntegrate, each point counts and so does their errors.

One way to avoid this is, as you already suggested, is to increase WorkingPrecision. Another way could be starting from a finite x, for example,

NIntegrate[BesselJ[9/2, x], {x, 10^-4, 1}]

0.000146199

Just for the record, at small x Integrate can also behave weired.

Integrate[BesselJ[9/2, x], {x, 0, 10^-4}] // N

-0.000048699

Integrate[BesselJ[9/2, x], {x, 0., 10^-4}] // N

0.

just by changing 0 to 0.