# Implementing Gauss–Laguerre quadrature for double integral of oscillatory function

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.

• @AntonAntonov, thank you, it work better now. I need to read up on those – Yuriy S Sep 8 '18 at 0:12
• I have also tried to implement Levin algorithm explicitly, but I have precision loss issues with LinearSolve, and the runtime is not very good – Yuriy S Sep 8 '18 at 1:23
• It's a good idea to include the error messages (I got one) that your code generates. In this case the one I get suggests that one result is somewhat inaccurate. – Michael E2 Sep 8 '18 at 2:17
• @MichaelE2, you are right, I didn't include it, but I mentioned the error. I'll make sure to include next time – Yuriy S Sep 8 '18 at 10:07
• I see. I thought you meant that the result was not consistent with a value known to you, not that NIntegrate has issued an error message. IMO, including the actual name NIntegrate::eincr potentially helps other users search for Q&As with similar problems. Also, I originally mentioned it because occasionally you get different messages or no messages in different versions of Mathematica, and I want to check whether you were getting what I was getting. – Michael E2 Sep 8 '18 at 12:56

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?

• Thank you very much! This seems to work very well for the double integral up to very high values of the parameters! I don't need to reinvent the wheel – Yuriy S Sep 8 '18 at 9:39
• @YuriyS Great, good luck! – Anton Antonov Sep 8 '18 at 9:52

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.

• Thank you. Turns out, I can change the variable to $s=r/b$ and then integrate only to $s \leq 1+\epsilon$, for example, $s \leq 5/4$, and the rest becomes a single integral with the square root approximated by its Taylor series. So I will also try your approach on the double integral with the limits $0\leq u \leq 1$ and $0 \leq s \leq 5/4$, see if it works faster – Yuriy S Sep 9 '18 at 14:13
• Turns out, your method performs better when the upper limit is $\infty$, but worse if the upper limit is $5/4$. In any case, your answer had been very helpful – Yuriy S Sep 9 '18 at 19:17