Fast Method for Numerically Integrating Bessel Functions over Finite Range

Question

I want to evaluate integrals of the kind

$$ \int_0^1 x \tanh{x}\, \mathcal{J}_{l + k}(\alpha_{n_1,l+k}\,x)\mathcal{J}_{l}(\alpha_{n_2,l}\,x) \, dx $$ where $\alpha_{n,l}$ is the $n^{th}$ zero of the $l^{th}$ Bessel function, $\mathcal{J}_l(x)$.

Currently, I'm using

Integ[l_, k_, n1_, n2_] := NIntegrate[x Tanh[x] BesselJ[l + k, x bjz[l + k, n1]] BesselJ[l, x bjz[l, n2]], {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}, WorkingPrecision -> 6, PrecisionGoal -> 5, AccuracyGoal -> 10];

where I define

bjz[n_?NumericQ, k_?NumericQ] := bjz[n, k] = N[BesselJZero[n, k]]

Let's say $k = 1$ here. For small n's and $l$, this method works fast enough for my purposes but say I take $l = 100, n_1 = 100, n_2 = 100$. Then it takes around 1 sec.

I tried using "LevinRule" and while that does speed up the integrals at large $l$ and $n$, it's slower than the other method at smaller values.

What would be the fastest way of evaluating such integrals such that it works for $0\leq l \leq 100$ and $0\leq n \leq 100$? I'm primarily interested in $k = 1,2$, although that shouldn't change much. There must be a general method for evaluating integrals of Bessel functions over a finite range, since these are pretty ubiquitous.

Answer 1

Well, there is no magical solution just because bessel-functions are involved.

But we can do the gold old brute-force-like testing. So lets define the measure Function:

MeasureTimeForInteg[l_?NumericQ,k_,n1_,n2_,method_]:=Timing[Integ[l,k,n1,n2,method];][[1]]

And because we don't like to interpret 3D we use a simple ListPlot approach:

MeasureTimingPlot[k_,n1_,n2_,maxL_:50]:=(
methods={Automatic,{Automatic,"SymbolicProcessing"->0},"LevinRule","GlobalAdaptive","ClenshawCurtisRule","GaussKronrodRule","LobattoKronrodRule"};
ListLinePlot[Transpose[ParallelTable[MeasureTimeForInteg[#,k,n1,n2,m],{m,methods}]&/@Range[1,maxL]],PlotRange->All,PlotLegends->(ToString/@methods),FrameLabel->{"l","time"},PlotLabel->{k,n1,n2}, PlotTheme -> "Scientific"]
)

You can already see, that i selected only a few methods which should be suitable in one or another way.

So we can look at the graphs for special values:

MeasureTimingPlot[1,1,1,70]
MeasureTimingPlot[1,100,1]
MeasureTimingPlot[1,20,20]
MeasureTimingPlot[1,50,50]

So we see, the best results in terms of time used comes from the "GlobalAdaptive"-Method while the plain Automatic and "LevinRule" also makes a decent job.

Therefore I would use GlobalAdaptive. It gives you the best performance out of them all.

Answer 2

The Gauss rule converges rapidly as the number of sample points increases, as does the Clenshaw-Curtis rule. I would try using "GaussBerntsenEspelidRule" with a heuristic that increases the number of points as the number of oscillations increases. This rule has a poor error estimator, so it's best to turn off recursion (MaxRecursion -> 0) and manage error manually. Below is a proof of concept:

integ[l_, k_, n1_, n2_] := 
  NIntegrate[
   x Tanh[x] BesselJ[l + k, x bjz[l + k, n1]] BesselJ[l, x bjz[l, n2]],
   {x, 0, 1}, MaxRecursion -> 0, 
   Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
     Method -> {"GaussBerntsenEspelidRule", 
       "Points" -> Round[9 + (l + k + n1 + n2)/2]}}];

Here is a test of the concept. Warning: It takes a long time to compute the "exact" integrals used for comparing the Gauss rule ones.

timings = ParallelTable[
   First@AbsoluteTiming[(Quiet@integ[l, k, n1, n2] - exact)/exact],
   {l, 1, 100, 33}, {k, 2}, {n1, 1, 100, 33}, {n2, 10, 100, 30}];
errors = ParallelTable[
   Block[{exact},
    exact = 
     NIntegrate[
      x Tanh[x] BesselJ[l + k, x BesselJZero[l + k, n1]] BesselJ[l, x BesselJZero[l, n2]],
      {x, 0, 1}, WorkingPrecision -> 32];
    (Quiet@integ[l, k, n1, n2] - exact)/exact
    ],
   {l, 1, 100, 33}, {k, 2}, {n1, 1, 100, 33}, {n2, 10, 100, 30}];

Here is a visualization of the timings and relative error:

GraphicsRow[{
  Histogram[res[[All, All, All, All, 1]] // Abs // Flatten, 
   PlotLabel -> "Timings"],
  Histogram[res[[All, All, All, All, 2]] // Abs // Flatten // Log10, 
   PlotLabel -> "Log10 rel. error"]}]

The maximum timing and relative error are:

timings // Max
(*  0.104274  *)

errors // Abs // Max
(*  6.94594*10^-9  *)