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

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.

• Since your $\epsilon_{mn}$ matrix is basically just the positive-frequency portion of the fast Fourier transform of the $\epsilon$ function, I feel like the use of 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. – DumpsterDoofus Aug 10 '14 at 20:31
• Question: Is your notation $\int_0^{2\pi}dp_xdp_y$ a shorthand for double integral $\int_0^{2\pi}\int_0^{2\pi}dp_xdp_y$ over a square region? – DumpsterDoofus Aug 10 '14 at 20:35
• @DumpsterDoofus; yes, that is a double integral, as you may have figured. Thanks. – ControleSam Aug 10 '14 at 21:29
• Yeah I think I can speed this up by a huge factor. I'll post in a while. – DumpsterDoofus Aug 10 '14 at 21:32
• Hmm, I found the reason why the FFT computation is so bad. According to aip.de/groups/soe/local/numres/bookcpdf/c13-9.pdf (from "Numerical Recipes in C"), it is tempting but very often massively inaccurate from a numerical perspective to compute Fourier coefficients using the FFT! Quoting from the book, "It is a sobering exercise to implement equation (13.9.6) for an integral that can be done analytically, and to see just how bad it is. We recommend that you try it." They then go on to provide a higher-order interpolation scheme. I'll give it a shot when I have time. – DumpsterDoofus Aug 11 '14 at 2:13

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...