0
$\begingroup$

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!

$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$
    – Michael E2
    Commented Apr 8, 2019 at 22:51
  • $\begingroup$ @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. $\endgroup$
    – Sid
    Commented Apr 9, 2019 at 7:27

1 Answer 1

2
$\begingroup$

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]

fig1

$\endgroup$
4
  • $\begingroup$ 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) $\endgroup$
    – Michael E2
    Commented Apr 9, 2019 at 22:07
  • $\begingroup$ @MichaelE2 I agree that it will be optimal PrecisionGoal ->3 and faster PrecisionGoal ->2. Contours vary slightly for each option. $\endgroup$ Commented Apr 9, 2019 at 23:33
  • $\begingroup$ 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. $\endgroup$
    – Sid
    Commented Apr 10, 2019 at 11:50
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Apr 10, 2019 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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