Why can't I change the value of MaxRecursion in NIntegrate when integrating BesselJ?

Question

Bug introduced in 8.0.4 or earlier and persisting through 11.0.1 or later

---

I am trying to evaluate this integral numerically $$ \int_0^{\infty } J_0(q R) \tanh(q) \, \mathrm{d}q $$ for large values of $R$. This makes the integrand oscillate more quickly and Mathematica gives incorrect answers. To deal with this I am trying to increase MaxRecursion in NIntegrate. Simply coding

 With[{R = 50},  
      NIntegrate[BesselJ[0, q R ] Tanh[q], {q, 0, ∞}, 
                 AccuracyGoal -> 12, PrecisionGoal -> 4, 
                 MaxRecursion -> 100]
 ]

throws no errors but it also does not increase computation time or give the correct answer.

If I set MinRecursion to a large value (larger than 9 - the default value in NIntegrate) in an attempt, I see an increase in computation time

 With[{R = 50},  
      NIntegrate[BesselJ[0, q R ] Tanh[q], {q, 0, ∞}, 
                 AccuracyGoal -> 12, PrecisionGoal -> 4,
                 MinRecursion -> 20, MaxRecursion -> 100]
 ]

I get an error saying

NIntegrate::minmax: MinRecursion (20) is greater than MaxRecursion (9).

I find this very confusing as I implicitly set the value of MaxRecursion in the code and it is not 9. Mathematica will allow my Min and Max Recursion if I delete the Bessel function and just have the Tanh in NIntegrate. My only thought is that this is some built-in property of BesselJ. Mathematica will also evaluate the BesselJ to arbitrary precision so I see no reason to limit the number of numerical subdivisions. Does anyone know a workaround?

P.S. Here is a some code which will quickly produce a plot of the integral as a function of $R$

 f[R_?NumericQ] := NIntegrate [BesselJ[0, q R ] Tanh[q], {q, 0, ∞}]

 LogLogPlot[
    f[R], {R, 1, 250}, PlotPoints -> 10, 
    MaxRecursion -> 1, AxesOrigin -> {0, 0}]

The code works up until $R$ is about 15 then gibberish for anything larger.

Answer 1

As has been noted by ruebenko in the comments, there does seem to be a bug in the handling of infinite-range Bessel function integrals when MinRecursion and MaxRecursion are both set to non-default values. For instance, even the simple

NIntegrate[BesselJ[0, x], {x, 0, ∞}, MinRecursion -> 10, MaxRecursion -> 15]

chokes with a NIntegrate::minmax error.

In any event, for the slightly more complicated

$$\int_0^\infty J_0(50u)\tanh\,u\,\mathrm du$$

what you can do is to explicitly use a method for infinite-range oscillatory integrals, and crank up WorkingPrecision while you're at it. For example, using Longman's method:

NIntegrate[BesselJ[0, 50 q] Tanh[q], {q, 0, ∞},
           Method -> "ExtrapolatingOscillatory", WorkingPrecision -> 90]
   2.1950746252821515546830074912679107125599945310570775933×10⁻³⁵

Hmm, a bit tiny. Is it actually zero? Let's check with something slightly different.

Let's take whuber's splitting suggestion. Using the identity

$$\tanh\,u=1-\exp(-u)\;\mathrm{sech}\,u$$

and exploiting the Hankel transform identity

$$\int_0^\infty J_0(cu)\,\mathrm du=\frac1{c},\quad c>0$$

we start by integrating the integral with $\mathrm{sech}$, using again Longman's method:

NIntegrate[BesselJ[0, 50 q] Exp[-q] Sech[q], {q, 0, ∞},
           Method -> "ExtrapolatingOscillatory", WorkingPrecision -> 90]
   0.01999999999999999999999999999999997804925374717848445316992508732089287121184172744576

which can be seen to be quite close to $1/50$. Subtracting this quantity from $0.02$, yields a result that agrees with the earlier attempt, so we now have a bit more trust in the results.

I had used Longman's method in these examples, but one could also have chosen to use the methods of Ooura-Mori ("DoubleExponentialOscillatory") or Levin ("LevinRule") instead.

Answer 2

Today I believe I encountered the same problem when trying to reproduce the result of this paper about Lamb's problem and solutions mentioned above doesn't help in my specific case. After struggling for a while I managed to find a work-around and I think it's worth sharing.

In short, if the integrate contains singular point(s) in addition, you may need to set higher WorkingPrecision and MaxRecursion, exclude the singular point(s) with Exclusions and explicitly set Method to "ExtrapolatingOscillatory". (Yeah, not one less.)

The following is the corresponding integral:

Clear[p, ξ, Nz, uztran]

(* the following is equation (16b) together with (11) etc. in the paper *)
Nz[p_] = With[{ν = 1/4, ρ = 7800, Ε = 210 10^9, r = 1, z = 0}, 
  With[{λ = (Ε ν)/((1 + ν) (1 - 2 ν)), μ = Ε/(2 (1 + ν))}, 
   With[{c = Piecewise[{{Sqrt[(λ + 2 μ)/ρ], # == 1}, {Sqrt[μ/ρ], # == 2}}] &}, 
    With[{k = p/c[#] &}, 
     With[{a = Sqrt[ξ^2 + k[#]^2]/k[2] &}, 
      With[{R = (2 ξ^2 + k[2]^2)^2 - 4 k[2]^2 a[1] a[2] ξ^2}, 
       (-(2 ξ^2 + k[2]^2) Exp[k[2] a[1] z] + 2 ξ^2 Exp[-k[2] a[2] z]) 
        (ξ k[2] a[1] BesselJ[0, ξ r])/R]]]]]];

(* the following is equation (15b) in the paper *)
uztran[p_, opt : OptionsPattern[]] := 
 (F gtran)/(2 π μ) NIntegrate[Rationalize[Nz[p], 0], {ξ, 0, ∞}, opt]

(* parameter that'll cause trouble, found by trial and error *)
specialp = 1.50026105063867726694999966600905294434039033506`12.86057158578721*^-10 - 
   1.67718225969505448435499738721675509494445398372`12.90898501143577*^-10 I;

uztran[specialp, 
  WorkingPrecision -> 64, 
  MaxRecursion -> 40, 
  Method -> "ExtrapolatingOscillatory", 
  Exclusions -> Denominator@Nz@specialp == 0] // AbsoluteTiming

{6.372237, 
  -(1/μ)(0.1193662073189147296405743133616902778179 + 
     7.5708093349255506613916610*10^-15 I) F gtran}

You can try modifying those 4 options to see what will happen.

P.S. It's worth emphasizing the necessity of Method -> "ExtrapolatingOscillatory" is actually another bug (this has been conformed by WRI): when you add Exclusions option to this integral, NIntegrate will spit out NIntegrate::bdmtd and NIntegrate::nsr without "ExtrapolatingOscillatory", at the moment it's the only work-around I found for this bug.