12
$\begingroup$

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]]

enter image description here

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 *)
$\endgroup$
2
  • 1
    $\begingroup$ NIntegrate[BesselJ[9/2, x], {x, 0, 1}, Method -> "LobattoKronrodRule"] gives better approximation $\endgroup$ Jun 8, 2016 at 8:53
  • 1
    $\begingroup$ Just as I noted in another thread, an expression in terms of SphericalBesselJ[] that Mathematica does not auto-expand to its trigonometric form will give better results: NIntegrate[Sqrt[2 x/π] SphericalBesselJ[4, x], {x, 0, 1}]. $\endgroup$ Jun 8, 2016 at 10:22

2 Answers 2

8
$\begingroup$

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.)

$\endgroup$
4
  • 3
    $\begingroup$ As an addendum: "SymbolicProcessing" -> 0 does absolutely nothing because BesselJ[] is already auto-expanded in the half-integer case, and the resulting expression is numerically iffy. $\endgroup$ Jun 8, 2016 at 12:04
  • $\begingroup$ @J.M. I should say I'm a little surprising, I didn't notice that NIntegrate evaluates its first argument while it owns HoldAll attribute, well, does NIntegrate always evaluate its first argument? (Maybe I should start a new question?) BTW @anton the // Expand is redundant :) $\endgroup$
    – xzczd
    Jun 9, 2016 at 3:33
  • $\begingroup$ I think, @xzczd, that the way to answer that would involve the IntegrationMonitor framework that Anton introduced to us previously. Which would be an interesting separate question, I would say. $\endgroup$ Jun 9, 2016 at 3:58
  • $\begingroup$ @j.m OK, I started a new question. $\endgroup$
    – xzczd
    Jun 9, 2016 at 6:09
2
$\begingroup$

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}]

enter image description here

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.

$\endgroup$

Your Answer

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

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