NIntegrate and Integrate of a Spherical Bessel function

Question

I am trying to integrate a spherical Bessel function. I have used both the Integrate and NIntegrate functions in Mathematica but the values given by each do not match. Any reason why this happens?

NIntegrate[((BesselJ[15/2, BesselJZero[15/2, 1]*r]))^2*r, {r, 0, a}]

NIntegrate::ncvb: "NIntegrate failed to converge to prescribed accuracy after 9 recursive
bisections in r near {r} = {3.13053*10^-59}. NIntegrate obtained
1.81126695652224073244637134881619178279603605269044825189376376942`65.954589770191*^2647
and 1.7853272819936654408129572054378095060361571646917142221312888537`65.954589770191*^2647
for the integral and error estimates."

1.811266956522241*10^2647

Integrate[((BesselJ[15/2, BesselJZero[15/2, 1]*r]))^2*r, {r, 0, a}]
0.

Edit: I forgot to mention $a=0.02$. I also noticed that the below integral should converge as I do not see any singularities but I'm not too familiar with spherical Bessel functions.

Integrate[((BesselJ[15/2, 0.0000001*r]))^2*r, {r, 0, a}]

Integrate::idiv: "Integral of
((-1.35135*10^47\r+1.7325*10^32\r^3-<<22>>\<<1>>+r^7)\<<1>>+<<1>>)^2/r^14 does not
converge on {0,0.02`}."

                   6                          -7
Integrate[6.3662 10  Power[Cos[12.5664 - 1. 10   r] - 

               47                    -7               32                    -7
     1.35135 10   Cos[12.5664 - 1. 10   r]   1.7325 10   Cos[12.5664 - 1. 10   r]
     ------------------------------------- + ------------------------------------ - 
                       6                                       4
                      r                                       r

            16                    -7                54                    -7
     3.78 10   Cos[12.5664 - 1. 10   r]   1.35135 10   Sin[12.5664 - 1. 10   r]
     ---------------------------------- - ------------------------------------- + 
                      2                                     7
                     r                                     r

             39                    -7             24                    -7
     6.237 10   Sin[12.5664 - 1. 10   r]   3.15 10   Sin[12.5664 - 1. 10   r]
     ----------------------------------- - ---------------------------------- + 
                      5                                     3
                     r                                     r

           8                    -7
     2.8 10  Sin[12.5664 - 1. 10   r]
     --------------------------------, 2], {r, 0, 0.02}]
                    r

Answer 1

Expressing the entire business in terms of SphericalBesselJ[] cures the problem:

NIntegrate[r With[{x = r BesselJZero[15/2, 1]}, 
                  Sqrt[2 x/π] SphericalBesselJ[7, x]]^2, {r, 0, 1/50}]
   1.1879560281974252*^-27

The nice thing about SphericalBesselJ[] is that it does not auto-evaluate to a potentially numerically unstable combination of trigonometric functions, unlike BesselJ[15/2, x].

Answer 2

It's just another story of precision. For the first sample:

a = 2/100;
NIntegrate[BesselJ[15/2, BesselJZero[15/2, 1] r]^2 r, {r, 0, a}, WorkingPrecision -> 40]
N[Integrate[BesselJ[15/2, BesselJZero[15/2, 1] r]^2 r, {r, 0, a}], 40]
(* Let's check the difference between the above two result *)
%% - %

1.187956028197538867114723184727333535080*10^-27
1.187956028197538867114723184727562859206*10^-27
-2.29324126*10^-58

For the second sample:

int2 = Integrate[BesselJ[15/2, 10^-7 r]^2 r, {r, 0, a}];
Block[{$MaxExtraPrecision = 290}, N[int2, 16]]

1.194601120645255*10^-148

To summarize, one should always keep the precision issue in mind when facing numerical calculation.

Answer 3

Try this:

    f = Interpolation[
  Table[{r, 10^12 BesselJ[15/2, BesselJZero[15/2, 1]*r]}, {r, 0, 0.02,
     0.0005}]];

and then integrate:

    10^-12*(NIntegrate[f[r], {r, 0, 0.02}] // Chop)

(*   1.67209*10^-14   *)

Have fun!

Later edit: to address the comment of xzczd: I, indeed, missed the square in my integration. If one adds this square and looks at the plot:

Plot[BesselJ[15/2, BesselJZero[15/2, 1]^2*r], {r, 0, 0.02}]

one can see the following:

It is easy to roughly estimate "by eye" the area under the curve as a fraction (by eye it is about 1/3) of the maximal height (0.00025), times the interval in which the curve is (visibly) different from zero (about 0.01):

 0.00025*0.01/3

(*  8.33333*10^-7   *)

yielding about 10^-6.

If one applies the approach I proposed above (with the function squared that I missed):

 f = Interpolation[
   Table[{r, 10^12 BesselJ[15/2, BesselJZero[15/2, 1]^2*r]}, {r, 0, 0.02, 0.0005}]];

10^-12*(NIntegrate[f[r], {r, 0, 0.02}] // Chop)

one gets:

(*   1.40079*10^-6   *)

in agreement with the above estimate, rather than 10^-27 or 10^-148.

Note that it is not necessary to multiply-divide by 10^12. This is simply one of the traces of my earlier trials to force Mathematica to integrate properly without tricks.

Finally, I think there is some problem with this integration, when it is being done without any trick. When one looks at the plot above, one sees that the function is smooth, non/oscillating, non-diverging in the interval of interest, and, therefore, there should be no problem to integrate it numerically directly. Mathematica should be able to take care about precision automatically, as it, indeed, does in other comparable cases. The integration MUST work without any trick. I guess that the bug (if any) is somewhere in the Bessel function, rather than in the integration itself.
The workaround, however, seems to work properly.

Have fun still!