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

Bug introduced in 8.0.4 or earlier and persists through 10.0.2.

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] := Module[{},  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.

Thanks.

-
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. – Szabolcs May 7 '13 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. – Szabolcs May 7 '13 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. – whuber May 7 '13 at 17:54
At least for your simpler problem, "ExtrapolatingOscillatory" (Longman's method) and "DoubleExponentialOscillatory" (Ooura-Mori method) both work well. – J. M. May 7 '13 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 '13 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.

-