Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
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
2  
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
3  
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
show 3 more comments

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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