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

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.

• Welcome to Mathematica.SE, and thank you for formatting the question for readability! Please do not add the [bugs] tag to new questions. It'll be added later if the consensus of the community is that this is indeed a bug. May 7, 2013 at 16:27
• I don't know what's happening here, but this page may be useful for you (if you're not yet aware of it). It details the various method options. May 7, 2013 at 16:33
• Depending on your accuracy needs, approximate $\tanh(z)\approx 1$ for $z$ sufficiently large, such as $z\gg 20$. This works because $1 - \tanh(z)$ decays exponentially. The integral of $J_0$ can be evaluated in closed form, so numerical integration is needed only for the product of the Bessel function with $\tanh(z)-1$ from $0$ to this small threshold. May 7, 2013 at 17:54
• At least for your simpler problem, "ExtrapolatingOscillatory" (Longman's method) and "DoubleExponentialOscillatory" (Ooura-Mori method) both work well. May 7, 2013 at 19:24
• I filed the issue that Max- and MinRecusion can not be set simultaneously as a bug. For the issue at hand (oscillatory integrand) you may want to try the "LevinRule". With[{R = 50}, NIntegrate[BesselJ[0, q R] Tanh[q], {q, 0, \[Infinity]}, AccuracyGoal -> 12, PrecisionGoal -> 4, Method -> "LevinRule"]] (*-1.72182*10^-15*). Another approach could be to increase the WorkingPrecision.
– user21
May 8, 2013 at 6:04

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.

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.5002610506386772669499996660090529443403903350612.86057158578721*^-10 -
1.6771822596950544843549973872167550949444539837212.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.

• I dunno... the integrand seems to be qualitatively different enough to warrant a separate thread about this, rather than have it as this answer. Apr 14, 2016 at 14:29
• @J.M. I did hesitate for a while but finally decided to post this answer. Anyway this integral contains BesselJ and if you tried something like uztran[specialp, WorkingPrecision -> 64, MinRecursion -> 20, MaxRecursion -> 40] you'll see the same warning as that mentioned in the question. Apr 14, 2016 at 14:37
• So many With[] constructions. Is it a good habit to use this?
– xyz
Apr 14, 2016 at 15:35
• @ShutaoTANG Mainly because I'm still in v9 and want to make the post self-contained so not that willing to use LetL etc. in this sample. You can have a look at this post: mathematica.stackexchange.com/a/54874/1871 Apr 15, 2016 at 2:28