7
$\begingroup$

"After multiplying the integrand of NIntegrate with -1, the Precision of the output will change." ← Sounds silly, huh? But this seems to be true at least for numerical integral internally using "ExtrapolatingOscillatory" method. Just try the following example:

Precision /@ 
 NIntegrate[{1, -1} BesselJ[0, x], {x, 0, ∞}, WorkingPrecision -> 32, 
  Method -> "ExtrapolatingOscillatory"]

{31.0265, 25.0279}

It's not necessary to set Method -> "ExtrapolatingOscillatory" manually in this sample, I added the option just to emphasize.

Of course in the above example the difference of precision is small and isn't a big deal, but in some cases the difference can be drastic, for example the following I encountered in this problem:

f[p_, ξ_] = -(5 p Sqrt[(5 p^2)/6 + ξ^2] )/(
  4 (-4 ξ^2 Sqrt[(5 p^2)/6 + ξ^2] Sqrt[(5 p^2)/2 + ξ^2] + ((5 p^2)/2 + 2 ξ^2)^2));

pmhankel[p_, sign_: 1, prec_: 32] := 
 NIntegrate[sign ξ BesselJ[0, ξ] f[p, ξ], {ξ, 0, ∞}, 
  WorkingPrecision -> prec, Method -> "ExtrapolatingOscillatory"]

preclst = Table[Precision@pmhankel[#, sign] & /@ Range@32, {sign, {1, -1}}]

ListLinePlot[preclst, PlotRange -> All]

enter image description here

It's not necessary to set Method -> "ExtrapolatingOscillatory" manually in this sample, I added the option just to emphasize.

How to understand the behavior? Except for calculating every integral twice and choosing the better one, how to circumvent the problem?

Improve this question
$\endgroup$
8
  • $\begingroup$ No. Doesn't sound silly to me. (1) NIntegrate does some symbolic pre-processing and anything you do you might affect that. That is reason enough to not find this surprising. If you haven't already, try ("SymbolicProcessing"->0). (2) And are you absolutely sure that changing the sign wouldn't significantly affect the evaluation for the numerical method you chose? That isn't obvious to me either. (reference.wolfram.com/language/tutorial/…) $\endgroup$
    – Searke
    Commented Apr 22, 2016 at 15:14
  • $\begingroup$ Yes, I wouldn't expect "ExtrapolatingOscillatory" to work without preliminary symbolic processing, how else can the zeroes where the integrand will be split at be determined? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Apr 23, 2016 at 6:06
  • $\begingroup$ @Searke Thanks for the response. (1) With "SymbolicProcessing->0" the problem remains: pmhankelTest[p_, sign_: 1, prec_: 16] := NIntegrate[sign ξ BesselJ[0, ξ] f[p, ξ], {ξ, 0, ∞}, WorkingPrecision -> prec, Method -> {"ExtrapolatingOscillatory", "SymbolicProcessing" -> 0}];pmhankelTest[32, #, 32] & /@ {-1, 1} (2) I tried IntegrationMonitor mentioned in this answer, the {"Boundaries", "Dimension", "Error", "GetRule", "Integrand"} etc. seems to be all the same, and the only difference between "Integral" is the sign. $\endgroup$
    – xzczd ♦
    Commented Apr 23, 2016 at 6:22
  • $\begingroup$ @J.M. I also tried manually implement "ExtrapolatingOscillatory", and the problem doesn't show up in my (much slower) implementation: zero[i_] := Piecewise[{{BesselJZero[0, i], i > 0}}];separatepmhankel[p_?NumericQ, sign : 1 | -1, i_?NumericQ, prec_] := NIntegrate[sign ξ BesselJ[0, ξ] f[p, ξ], {ξ, zero@i, zero[i + 1]}, WorkingPrecision -> prec, MaxRecursion -> 40]; manualpmhankel[p_, sign_: 1, prec_: 16] := NSum[separatepmhankel[p, sign, i, prec], {i, 0, Infinity}, Method -> "AlternatingSigns", WorkingPrecision -> prec]; manualpmhankel[32, #, 32]&/@{1,-1} // AbsoluteTiming $\endgroup$
    – xzczd ♦
    Commented Apr 23, 2016 at 6:27
  • 1
    $\begingroup$ @Xavier Interesting. Using the code in the update of this answer, I found that in this case NIntegrate internally switches to "LevinRule", and the output is indeed the same as that with option , Method -> "LevinRule". How about giving an answer? $\endgroup$
    – xzczd ♦
    Commented Apr 25, 2016 at 5:54

2 Answers 2

Reset to default
6
$\begingroup$

The oddity in this case comes from NSum which is being called in a certain way from NIntegrate. This is a simple example that has roughly the same behavior (note in this case the exact result is known to be $\mp \ln 2$):

NSum[(-1)^n/n, {n, 1, Infinity}, 
          Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, WorkingPrecision -> 32]

(* -0.6931471805599453094172318803247 *)

NSum[-(-1)^n/n, {n, 1, Infinity}, 
          Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, WorkingPrecision -> 32]

(* 0.693147180559945309417232 *)

where the second result has several digits fewer than the first.

Is that a bug? Not necessarily, because both results have at least 16 correct digits which certainly attains the default PrecisionGoal, which is WorkingPrecision/2.

Still, I agree the consistency could be improved in this case and I have filed a report for the developers to take a look.

Update

This has been improved in the just released Mathematica 11.0.

Improve this answer
$\endgroup$
7
$\begingroup$

Short answer:

One workaround is to use Method -> "LevinRule" instead.

Long answer:

As mentioned by Xavier in a comment above, changing BesselJ[0, x] to BesselJ[0, Re[x]] resolves the issue:

NIntegrate[{1, -1} BesselJ[0, Re@x], {x, 0, ∞}, WorkingPrecision -> 32, 
 Method -> "ExtrapolatingOscillatory"]
    Precision /@ %

{0.99999999999999999999999999999979, -0.99999999999999999999999999999979}
{32., 32.}

But why it works? Codes in the UPDATE of this answer solve the mystery. The truth is, NIntegrate has internally switched to "LevinRule":

enter image description here

Pictured by Simon Wood's shadow.

The output is indeed the same as that with option Method -> "LevinRule:

NIntegrate[{1, -1} BesselJ[0, x], {x, 0, ∞}, WorkingPrecision -> 32, 
 Method -> "LevinRule"]
Precision /@ %

{0.99999999999999999999999999999979, -0.99999999999999999999999999999979}
{32., 32.}

Notice in this case you'll be unable to adjust the MaxRecursion option because of a bug mentioned here, but this seems not to be a big deal:

pmhankelLevin[p_, sign_: 1, prec_: 32] := 
 NIntegrate[sign ξ BesselJ[0, ξ] f[p, ξ], {ξ, 0, ∞}, 
  WorkingPrecision -> prec, Method -> "LevinRule"]

lst1 = pmhankel[#, -1] & /@ Range@32;
(* You'll see warning NIntegrate::ncvb and be unable to eliminate it 
   by adjusting MaxRecursion option because of the bug mentioned above 
   when generating lst2 and lst3 *)
lst2 = pmhankelLevin /@ Range@32;
lst3 = pmhankelLevin[#, -1] & /@ Range@32;
(* But the difference between lst2, lst3 and lst1 is negligible: *)
Max /@ {(lst1 + lst2)/lst1, (lst1 - lst3)/lst1}

{2.27320119973*10^-6, 2.27320119973*10^-6}

Though a workaround is found, the reason why the precision is influenced by the sign when "ExtrapolatingOscillatory" method is used still remains unclear. I'm looking forward to answer(s) addressing the issue.

Improve this answer
$\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.