How can I reduce computation time while still obtaining a good approximation for my function?

Question

I am new to any CAS (and Mathematica, for that matter) and new to StackExchange too, so forgive me and correct me on any mistakes.

I have this function: $J_p=\sum_{m,n=1}^{\infty} \epsilon_{mn}f_{mn}\sum_{k=-\infty}^{\infty}\frac{J_k^2(\beta)(m\Omega+k\omega)}{1+(m\Omega+k\omega)^2}$ where $\epsilon_{mn}=-\frac{m n}{4\pi^2}\int_0^{2\pi}\epsilon(p_x,p_y)\exp(-i(m p_x+n p_y))\,dp_x dp_y$ and $f_{mn}=-\frac{m n}{4\pi^2}\int_0^{2\pi}\frac{\exp(-i(m p_x+n p_y))}{1+\exp(-\epsilon(p_x,p_y))}\,dp_x dp_y$ where again $\epsilon(p_x,p_y)=\sqrt{1+4\cos\left(\frac{p_y}{2}\right)\cos\left(\frac{p_x\sqrt{3}}{2}\right)+4\cos^2\left(\frac{p_y}{2}\right)}$.

Here is my Mathematica code to evaluate this:

Off[NIntegrate::ncvi];
epsilonCoeffsMMA[cl_] := Module[{reComp, imComp},
   reComp[m_, n_] := (-m n)/(4 \[Pi]^2)
      NIntegrate[
      Re[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)
         Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]}, 
      Method -> "Trapezoidal", MaxRecursion -> 100];
   imComp[m_, n_] := (-m n)/(4 \[Pi]^2)
      NIntegrate[
      Im[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)
         Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]}, 
      Method -> "Trapezoidal", MaxRecursion -> 100];
   emnMatrix = Table[0, {m, 1, cl}, {n, 1, cl}];
   Do[emnMatrix[[m, n]] = reComp[m, n] + I imComp[m, n], {m, 1, 
     cl}, {n, 1, cl}];
   ];
boltzECoeffsMMA[cl_] := Module[{reComp, imComp},
reComp[m_, n_] := (-m n)/(4 \[Pi]^2)
      NIntegrate[
      Re[Exp[-I (m px + n py)]/(1 + 
          Exp[-(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 
              4 Cos[py/2]^2)^(1/2)])], {px, 0, 2 \[Pi]}, {py, 0, 
       2 \[Pi]}, Method -> "Trapezoidal", MaxRecursion -> 100];
   imComp[m_, n_] := (-m n)/(4 \[Pi]^2)
      NIntegrate[
      Im[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)
         Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]}, 
      Method -> "Trapezoidal", MaxRecursion -> 100];
   fmnMatrix = Table[0, {m, 1, cl}, {n, 1, cl}];
   Do[fmnMatrix[[m, n]] = reComp[m, n] + I imComp[m, n], {m, 1, 
     cl}, {n, 1, cl}];
   ];
jPMMA[coeffLim_, kernLim_] := Module[
  {cl = coeffLim, kl = kernLim, px, py},
  epsilonCoeffsMMA[cl]; boltzECoeffsMMA[cl];
  coeffMatrix = emnMatrix fmnMatrix;
  sumMatrix = 
   Table[Sum[(
     BesselJ[k, \[Beta]]^2 (m \[CapitalOmega] + k \[Omega]))/(
     1 + (m \[CapitalOmega] + k \[Omega])^2), {k, -kl, kl}], {m, 1, 
     cl}, {n, 1, cl}];
  jParaMMA = Total[coeffMatrix sumMatrix, 2];
  ];

This generates a function jParaMMA which I can Plot after I have made the call jPMMA[a,b] for some integers; a and b. For example

jPMMA[10, 10];
Plot[Evaluate@
Re[jParaMMA /. {\[Beta] ->2, \[Omega] -> {0, 2, 4, 6, 8}}], {\[CapitalOmega], 0, 20}, PlotRange -> Full]

for which

First[Timing[jPMMA[10, 10]]]

gives

115.437500

My question is: How can I obtain similar results, possibly with more terms (i.e. from running jPMMA[50, 60], say.) in a shorter time,? Thank you.

PS: I used the Off[NIntegrate::ncvi] because I do not know how to eliminate it from my numerical integration and I'd be glad to obtain some help for that too. Also, I used the Trapezoidal method because I noticed it gave a faster approximation even when coupled with MaxRecursion -> 100. I have tried with the Cuba library implementation in both Mathematica and Maple, which I was led to by this post, and the approximations are appreciably close.

Answer 1

Timing under 20 seconds on my computer now. Ok, your original program took about 60 seconds on my computer meaning that my computer is faster.

The dramatical gain of time is due to halfing the MaxRecursion option value. The plot still shows no visible difference.

I replaced Pi-Symbol by Pi for increasing readability in forum. I tested some scenarios, and replaced your initialization with 0 for the matrix and made px, py local, here you can't win seconds, but shortens progam. Here I tested also with ClearSystemCache and AbsoluteTiming

your Timing result on my computer as image

After the changes mentioned the result is:

Here my changes:

Clear@"`*";
Off[NIntegrate::ncvi];
maxRecursion = 50;
epsilonCoeffsMMA[cl_] := 
  Module[{reComp, imComp, px, py}, 
   reComp[m_, 
     n_] := (-m n)/(4 Pi^2) NIntegrate[
      Re[(1 + 4 Cos[py/2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/
           2) Exp[-I (m px + n py)]], {px, 0, 2 Pi}, {py, 0, 2 Pi}, 
      Method -> "Trapezoidal", MaxRecursion -> maxRecursion];
   imComp[m_, 
     n_] := (-m n)/(4 Pi^2) NIntegrate[
      Im[(1 + 4 Cos[py/2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/
           2) Exp[-I (m *px + n*py)]], {px, 0, 2 Pi}, {py, 0, 2 Pi}, 
      Method -> "Trapezoidal", MaxRecursion -> maxRecursion];
   emnMatrix = 
    Table[reComp[m, n] + I imComp[m, n], {m, 1, cl}, {n, 1, cl}]];
boltzECoeffsMMA[cl_] := 
  Module[{reComp, imComp, px, py}, 
   reComp[m_, 
     n_] := (-m n)/(4 Pi^2) NIntegrate[
      Re[Exp[-I (m *px + n * py)]/(1 + 
          Exp[-(1 + 4 Cos[py/2] Cos[(px Sqrt[3])/2] + 
               4 Cos[py/2]^2)^(1/2)])], {px, 0, 2 Pi}, {py, 0, 2 Pi}, 
      Method -> "Trapezoidal", MaxRecursion -> maxRecursion];
   imComp[m_, 
     n_] := (-m n)/(4 Pi^2) NIntegrate[
      Im[(1 + 4 Cos[py/2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/
           2) Exp[-I (m*px + n*py)]], {px, 0, 2 Pi}, {py, 0, 2 Pi}, 
      Method -> "Trapezoidal", MaxRecursion -> maxRecursion];
   fmnMatrix = 
    Table[reComp[m, n] + I imComp[m, n], {m, 1, cl}, {n, 1, cl}]];
jPMMA[coeffLim_, kernLim_] := 
  Module[{cl = coeffLim, kl = kernLim}, epsilonCoeffsMMA[cl]; 
   boltzECoeffsMMA[cl];
   coeffMatrix = emnMatrix * fmnMatrix;
   sumMatrix = 
    Table[Sum[(BesselJ[k, β]^2 (m Ω + 
           k ω))/(1 + (m Ω + 
            k ω)^2), {k, -kl, kl}], {m, 1, cl}, {n, 1, cl}];
   jParaMMA = Total[coeffMatrix * sumMatrix, 2];];

So it's up to you, if tuning MaxRecursion is worth a compromise. Pick what you like, just some ideas...