Fast Method for Numerically Integrating Bessel Functions over Finite Range

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]

enter image description here enter image description here enter image description here enter image description here

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.


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"]}]

Mathematica graphics

The maximum timing and relative error are:

timings // Max
(*  0.104274  *)

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

Tags:

Numerical Integration

Performance Tuning

Special Functions

Recent Posts