Implementing Gauss–Laguerre quadrature for double integral of oscillatory function

Question

I need an efficient way to evaluate a certain double integral for a wide range of parameters:

$$I=(\gamma_l \gamma_L )^{3/2} \int_0^\infty \int_0^{1} \frac{r ~d r ~dx}{\sqrt{r^2+b^2 x^2}} \times \\ \times L_l \left( 2 \gamma_l r \right) L_L \left(2 \gamma_L r \right) e^{-(\gamma_l+\gamma_L) r} \left((-1)^{(j-J)/2} \cos \beta^{-}_{jJ} x-(-1)^{(j+J)/2} \cos \beta^{+}_{jJ} x\right)$$

Where $j,J,l,L$ are integers, such that $j \pm J$ is even, and: $$\beta^{\pm}_{jJ}=\frac{\pi (j \pm J)}{2} \\ \gamma_l=\frac{1}{l+1/2}$$

while $b>1$ is just a real parameter. I may need to have $j,J$ as large as $100$ at least. $l,L$ can be under $20$.

A special case $j=J$ allows one to take half of the integral over $x$ exactly, however the second part still needs to be evaluated numerically.

I have tried to make a change of variables and evaluate function of the form:

$$F(\beta,s)=s \int_0^1 \frac{\cos \beta u~du}{\sqrt{s^2+u^2}}$$

However, that wasn't successful, so I decided to use Gauss-Laguerre quadrature for the integral over $r$, because it seems to fit well. I have obtained the formula:

$$I \approx \frac{(\gamma_l \gamma_L )^{3/2}}{(\gamma_l+\gamma_L)(N+1)^2} \sum_{n=1}^N \frac{L_l \left( \frac{2 \gamma_l}{\gamma_l+\gamma_L} r_n \right) L_L \left(\frac{2 \gamma_L}{\gamma_l+\gamma_L} r_n \right)~r_n}{L_{N+1}^2(r_n)} \times \\ \times \int_0^{1} \frac{(-1)^{(j-J)/2} \cos \beta^{-}_{jJ} x-(-1)^{(j+J)/2} \cos \beta^{+}_{jJ} x}{\sqrt{1+\frac{b^2}{(\gamma_l+\gamma_L)^2 r_n^2} x^2}} dx$$

This seems to work for small $j,J$, but both performance and accuracy suffer when compared to Mathematica double integration routine. For example (updated after a comment about using PrecisionGoal):

j = 11;
J = 23;
l = 2;
L = 5;
b = 10;
a1 = N[1/(l + 1/2), 20];
a2 = N[1/(L + 1/2), 20];
b1 = N[(\[Pi] j)/2, 20];
b2 = N[(\[Pi] J)/2, 20];
bm = b1 - b2;
bp = b1 + b2;
j1 = If[EvenQ[(j - J)/2], 1, -1];
j2 = If[EvenQ[J], 1, -1];
A0 = j1  N[(a1 a2)^(3/2), 20];
Nm = 10 Max[l, L, j, J];
A = (j1 (a1 a2)^(3/2))/((a1 + a2) (Nm + 1)^2);
S0 = (a1 + a2)^2 b^2;
S1 = b^2;
R = x /. NSolve[LaguerreL[Nm, x] == 0, 20];

T1 = AbsoluteTiming[R1 = (2 a1)/(a1 + a2) R;
  R2 = (2 a2)/(a1 + a2) R;
  w = (A R)/LaguerreL[Nm + 1, R]^2 LaguerreL[l, R1] LaguerreL[L, R2];
  S = S0/(R^2) ;
  If[EvenQ[j + J],   
   NIntegrate[
    Total[ w/Sqrt[1 + S u^2]] (Cos[bm u] - j2 Cos[bp u]), {u, 0, 1}, 
    PrecisionGoal -> 5, Method -> "LevinRule"], 0]]
T2 = AbsoluteTiming[
  If[EvenQ[j + J], 
   A0 NIntegrate[(
     LaguerreL[l, 2 a1 r] LaguerreL[L, 
       2 a2 r] Exp[-r (a1 + a2)] r (Cos[bm u] - j2 Cos[bp u]))/Sqrt[
     r^2 + S1 u^2], {u, 0, 1}, {r, 0, Infinity}, PrecisionGoal -> 5], 
   0]]

I get:

{0.295737, 0.000107986}
{26.3857, 0.000107935}

For the double integral there's also extreme increase of error for GlobalAdaptive method (NIntegrate::eincr). Monte Carlo works better, but still takes very long.

---

As you can see, quadrature rule is more computationally stable and much faster*, but doesn't give very accurate results for a small number of nodes. Meanwhile, double integration suffers from both singularity and oscillations, which leads to loss of precision.

For larger $j,J$ double integration is extremely slow, and I'm not sure if it even gives results. However, quadrature is slow as well and only converges for large number of nodes.

My questions are: is there a way to improve my implementation or find some other way of computing this integral for the specified range of parameters, with some degree of accuracy?

Additionally, if I still use the quadrature rule, how do I choose the number of nodes correctly?

In any case, for large $j,J$ double integral doesn't evaluate in several minutes and that's unacceptable. However, how do I check the accuracy of the approximate method without knowing the exact value?

I'm sorry for the length of the post, thank you for your time!

---

* that if we don't include the evaluation of the roots for the Laguerre polynomial. The roots of course can be pre-computed.

Answer 1

It seems it is faster to use Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}:

T1 = AbsoluteTiming[R1 = (2 a1)/(a1 + a2) R;
  R2 = (2 a2)/(a1 + a2) R;
  w = (A R)/LaguerreL[Nm + 1, R]^2 LaguerreL[l, R1] LaguerreL[L, R2];
  S = S0/(R^2);
  If[EvenQ[j + J], 
   NIntegrate[
    Total[w/Sqrt[1 + S u^2]] (Cos[bm u] - j2 Cos[bp u]), {u, 0, 1}, 
    PrecisionGoal -> 8, 
    Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", 
      "SymbolicProcessing" -> 0}], 0]]

(* {0.221031, 0.000107986} *)


T2 = 
 AbsoluteTiming[
  If[EvenQ[j + J], 
   A0 NIntegrate[(LaguerreL[l, 2 a1 r] LaguerreL[L, 
         2 a2 r] Exp[-r (a1 + a2)] r (Cos[bm u] - j2 Cos[bp u]))/
      Sqrt[r^2 + S1 u^2], {u, 0, 1}, {r, 0, Infinity}, 
     PrecisionGoal -> 8, 
     Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", 
       "SymbolicProcessing" -> 0}], 0]]

(* {1.05594, 0.000107929} *)

Are the results something you expect?

Answer 2

This is a comment that extends Anton's answer.

Outside any small neighborhood of the origin, the integrand is analytic and you can take advantage of the spectral convergence of the Gauss and Gauss-Kronrod rules by increasing the number of Gauss points, especially along the u-axis. Despite the r-integration being out to infinity, the convergence rate seems to be aided by the exponential decay of the integrand. For the u-integration, convergence should begin around "Points" -> 0.5 bp or a little higher, and it reaches the desired goal a little under bp points. The r-integration is a little harder to predict (for me). A modestly high setting such as "Points" -> 21 greatly speeds up the calculation, and the optimal setting for the given example is "Points" -> 51.

T1 = AbsoluteTiming[
  R1 = (2 a1)/(a1 + a2) R;
  R2 = (2 a2)/(a1 + a2) R;
  w = (A R)/LaguerreL[Nm + 1, R]^2 LaguerreL[l, R1] LaguerreL[L, R2];
  S = S0/(R^2);
  If[EvenQ[j + J],
   NIntegrate[
    Total[w/Sqrt[1 + S u^2]] (Cos[bm u] - j2 Cos[bp u]), {u, 0, 1}, 
    PrecisionGoal -> 8, 
    Method -> {"GaussKronrodRule", "Points" -> Ceiling[0.8 bp]}], 0]]
(*  {0.095939, 0.000107986}  *)

T2 = AbsoluteTiming[
  If[EvenQ[j + J], 
   A0 NIntegrate[
     (LaguerreL[l, 2 a1 r] LaguerreL[L, 2 a2 r] Exp[-r (a1 + a2)] r *
      (Cos[bm u] - j2 Cos[bp u]))/Sqrt[r^2 + S1 u^2],
     {u, 0, 1}, {r, 0, Infinity}, 
     PrecisionGoal -> 5, 
     Method -> {
       {"GaussKronrodRule", "Points" -> Round[0.8 bp]}, 
       {"GaussKronrodRule", "Points" -> 51}}], 0]]
(*  {0.033414, 0.000107929}  *)

Another approach that is a little worse but still a significant improvement is to use the options:

Method -> {"GaussKronrodRule", "Points" -> 21},
MinRecursion -> Floor[Log2[bp/21]]

The number 21 can be replace by a lower or higher number. Forcing preliminary recursive subdivision with MinRecursion speeds up the default setting 5 quite a bit as well. Ceiling sometimes to yield a faster integration than Floor and vice versa.