6
$\begingroup$

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.

$\endgroup$
7
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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