Parallelise numerical integral for contour plot

I have the following integral that I would like to evaluate:

\[Tau]=0.9; \[Sigma]=0.9; a=13; b=8; c=11; R1=1.2; I1=0.8; R2=3.4; I2=3;
A[x_, y_] := 2/\[Pi] ((-1)^x y^x)/(2 - y)^(1 + x);

Pintegrand[m_, n_, r1_, r2_, \[Theta]1_, \[Theta]2_] := (A[m, \[Tau]] A[n, \[Sigma]])/(\[Pi]^2) r1 r2 LaguerreL[m, (4 r1^2)/(2 - \[Tau])] LaguerreL[n, (4 r2^2)/(2 - \[Sigma])] Exp[-((2 \[Tau] r1^2)/(2 - \[Tau]))] Exp[-((2 \[Sigma] r2^2)/(2 - \[Sigma]))] Exp[- (1/2) (-R2 + r2 Cos[\[Theta]2])^2 + (-I1 + r1 Sin[\[Theta]1]) (-((b (-R1 + r1 Cos[\[Theta]1]))/(-b^2 + a c)) + (a (-I1 + r1 Sin[\[Theta]1]))/(-b^2 + a c)) + (-R1 + r1 Cos[\[Theta]1]) ((c (-R1 + r1 Cos[\[Theta]1]))/(-b^2 + a c) - (b (-I1 + r1 Sin[\[Theta]1]))/(-b^2 + a c)) + (-I2 + r2 Sin[\[Theta]2])^2] // Simplify;

P11 = Integrate[Pintegrand[1, 1, r1, r2, \[Theta]1, \[Theta]2], {r1, 0, \[Infinity]}, {\[Theta]1, 0, 2 \[Pi]}, {r2, 0, \[Infinity]}, {\[Theta]2, 0, 2 \[Pi]}]

The integrand function is the integrand of the integral, which is a Gaussian function weighted with the Laguerre polynomials (which is what is making the calculation lengthy). The variables r1, \[Theta]1, r2 and \[Theta]2 are integration variables.

I would like to generate a list of results with different values of $$m$$ and $$n$$ - that is P01, P10, P11, P12, P21, P22, P13, etc (first number corresponds to value of $$m$$ and second to $$n$$), with both $$m$$ and $$n$$ integers that each go from 0 to 10. I would like to use the result to create a contour plot of the result such as a ListContourPlot.

There are two post that request a similar objective, such as this and this. However the solutions appear to be specific to the integral in question and make use of different strategies, such as Parallelize, ParallelTable, and ParallelMap.

How can I most efficiently parallelise the numerical integral to generate a contour plot? Any help is appreciated!

• There are many undefined parameteris. It makes it difficult to test out ideas. Can you give a typical example of the kind of contour plot you wish to compute? – Michael E2 Apr 8 at 22:51
• @MichaelE2, thanks for your comment! I have defined the free parameters that are realistic values for the work. The remaining parameters are integration parameters. I've also explained the type of contour plot I'd like. – Sid Apr 9 at 7:27

The study of the function shows that here it is possible to limit it to the parameters specified in the intervals 0<=r1<=9,0<=r2<=9,0<=m<=10,0<=n<=10. To display a list of integrals, high accuracy is not needed, but it is desirable to use a logarithmic scale.

\[Tau] = 9/10; \[Sigma] = 9/10; a = 13; b = 8; c = 11; R1 =
12/10; I1 = 8/10; R2 = 34/10; I2 = 3;
A[x_, y_] := 2/\[Pi] ((-1)^x y^x)/(2 - y)^(1 + x);
Pintegrand[m_, n_, r1_,
r2_, \[Theta]1_, \[Theta]2_] := (A[m, \[Tau]] A[
n, \[Sigma]])/(\[Pi]^2) r1 r2 LaguerreL[
m, (4 r1^2)/(2 - \[Tau])] LaguerreL[
n, (4 r2^2)/(2 - \[Sigma])] Exp[-((2 \[Tau] r1^2)/(2 - \[Tau]))] \
Exp[-((2 \[Sigma] r2^2)/(2 - \[Sigma]))] Exp[-(1/2) (-R2 +
r2 Cos[\[Theta]2])^2 + (-I1 +
r1 Sin[\[Theta]1]) (-((b (-R1 + r1 Cos[\[Theta]1]))/(-b^2 +
a c)) + (a (-I1 + r1 Sin[\[Theta]1]))/(-b^2 +
a c)) + (-R1 +
r1 Cos[\[Theta]1]) ((c (-R1 + r1 Cos[\[Theta]1]))/(-b^2 +
a c) - (b (-I1 + r1 Sin[\[Theta]1]))/(-b^2 + a c)) + (-I2 +
r2 Sin[\[Theta]2])^2];

lst =
ParallelTable[{m, n,
Log[Quiet[
NIntegrate[
Pintegrand[m, n, r1, r2, \[Theta]1, \[Theta]2], {r1, 0,
9}, {\[Theta]1, 0, 2 \[Pi]}, {r2, 0, 9}, {\[Theta]2, 0,
2 \[Pi]}, WorkingPrecision -> 4]]]}, {m, 0, 10}, {n, 0,
10}]; // AbsoluteTiming

ListContourPlot[Flatten[lst, 1], Contours -> 20,
PlotRange -> All, ColorFunction -> Hue, FrameLabel -> {"m", "n"},
PlotLegends -> Automatic] • Using machine precision and PrecisionGoal -> 2 or 3 is faster than WorkingPrecision -> 4m which sets the PrecisionGoal to 2 (probably) and uses slower floating-point libraries to do the arbitrary-precision calculations. (+1) – Michael E2 Apr 9 at 22:07
• @MichaelE2 I agree that it will be optimal PrecisionGoal ->3 and faster PrecisionGoal ->2. Contours vary slightly for each option. – Alex Trounev Apr 9 at 23:33
• The r1 and r2 integration variables are really the radius in polar coordinates. Is the function not well behaved above 9? Also, could you explain the use of Log and Quiet in the ParallelTable. – Sid Apr 10 at 11:50
• @Sid When r1>9, r2>9, the function Pintegrand is close to zero. In this case, we can trim the integral. In case messages may appear, we use Quiet to clip messages. The logarithm is used to display at one scale, since the integral is changed by 12 orders of magnitude. – Alex Trounev Apr 10 at 12:50