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.
NIntegrate
might actually not even be necessary, and a much faster method may be available. I'll see what I can do and post an answer in a few hours if I can figure it out. $\endgroup$