Systematically Optimizing Code for Better Runtime

Question

Can someone please tell me how can I learn to systematically optimize my Mathematica code to be faster? I have seen a blog about the same but have no idea where to begin and how to implement this.

For example :- The following code takes more than 24 hours to run (actually it's still running, expected 30+ hours)

Ef[a_] := Pi^2*(a + 2)^2; 
Eb[a_] := Pi^2*(a + 1)^2; 
f[n_, x_] := Sqrt[2/((n + 2)^2 - 1)]*((n + 2)*Cos[Pi*(n + 2)*x] - Cot[Pi*x]*Sin[(n + 2)*Pi*x]); 
b[n_, x_] := Sqrt[2]*Sin[(n + 1)*Pi*x]; 
xf[m_, n_] := If[Mod[m - n, 2] == 0, 0, Integrate[f[n, x]*f[m, x]*x, {x, 0, 1}, 
     Assumptions -> {Element[n, Integers], Element[m, Integers]}]]; 
xb[m_, n_] := Integrate[b[n, x]*b[m, x]*x, {x, 0, 1}, Assumptions -> {Element[n, Integers], Element[m, Integers]}]; 
xt[m_, n_] := If[m == 0 || n == 0, xb[m, n], (1/2)*(xb[m, n] + xf[m - 1, n - 1])]; 
Z[T_] := Sum[E^(-(Eb[i]/T)), {i, 0, 10}]; 
y1[m_, t_] := -Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]; 
Y1[T_, t_] := (-(1/Z[T]))*Sum[Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]/E^(Eb[m]/T), 
     {m, 0, 10}]; 
file = OpenAppend["susypotwell_Y1T0.1.dat"]; 
Table[Export[file, {{t, Y1[0.1, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]
file = OpenAppend["susypotwell_Y1T1.dat"]; 
Table[Export[file, {{t, Y1[1, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]
file = OpenAppend["susypotwell_Y1T10.dat"]; 
Table[Export[file, {{t, Y1[10, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]

One way to do this could be to evaluate Y1[T,t] and define a new function using the output which would then be used in the Export expressions instead of Y1. Please note that here I am exporting data to plot because Mathematica does not automatically save data of plots if it needs to be modified later on. If Plot can somehow do this faster then I have no problem with it as well and I will just use this or similar to save the plot data inside notebook for later manipulations.

Furthermore, if say, one has optimized their code as much as it can be and it is still taking 24+ hours or something like that then what are the options they have?

Apologies for such a broad question. Even quick tips and suggestions would help a lot.

Edit :-

Step-1 : Use #-& notation for pure functions. (source)

Answer 1

Instead of 30 hours, do the job in 2.6 seconds.

Edit Used NIntegrate, as @flinty recommended.

    ClearAll["Global`*"]
(Ef[a_] = Pi^2*(a + 2)^2;
Eb[a_] = Pi^2*(a + 1)^2;
f[n_, x_] = 
  Sqrt[2/((n + 2)^2 - 1)]*((n + 2)*Cos[Pi*(n + 2)*x] - 
  Cot[Pi*x]*Sin[(n + 2)*Pi*x]);
b[n_, x_] = Sqrt[2]*Sin[(n + 1)*Pi*x];

(*    Table[xf[m, n] = 
  If[Mod[m - n, 2] == 0, 0, 
 Integrate[f[n, x]*f[m, x]*x, {x, 0, 1}]], {m, 0, 10}, {n, 0, 10}];
Table[xb[m, n] = Integrate[b[n, x]*b[m, x]*x, {x, 0, 1}], {m, 0, 
10}, {n, 0, 10}];   *)

Table[xf[m, n] = 
  If[Mod[m - n, 2] == 0, 0, 
NIntegrate[f[n, x]*f[m, x]*x, {x, 0, 1}, MaxRecursion -> 50]], {m,
0, 10}, {n, 0, 10}];
Table[xb[m, n] = 
  If[(1/2 (2 + m + n)) \[Element] Integers && m != n, 0, 
NIntegrate[b[n, x]*b[m, x]*x, {x, 0, 1}]], {m, 0, 10}, {n, 0, 
10}];

xt[m_, n_] = 
  If[m == 0 || n == 0, xb[m, n], (1/2)*(xb[m, n] + xf[m - 1, n - 1])];
Z[T_] = Sum[E^(-(Eb[i]/T)), {i, 0, 10}] // Simplify;
y1[m_, t_] = -Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*
  Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]; 
Y1[T_, t_] := (-(1/Z[T]))*
  Sum[Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*
    Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]/E^(Eb[m]/T), {m, 0, 10}];
tab1 = Table[{t, Y1[0.1, t]}, {t, -1., 1., 0.01}];
tab2 = Table[{t, Y1[1, t]}, {t, -1., 1., 0.01}];

tab3 = Table[{t, Y1[10, t]}, {t, -1., 1., 0.01}];
{ListLinePlot[tab1, Epilog -> {Red, Point@tab1}], 
 ListLinePlot[tab2, Epilog -> {Red, Point@tab2}], 
 ListLinePlot[tab3, Epilog -> {Red, Point@tab3}]}
) // Timing

Let me say it repeadetly: Avoid SetDelayed (:=) whereever you can. My opinion.

Answer 2

In order to determine where to concentrate your efforts, you need to know where your bottlenecks are.

To do this work through a single Y1 calculation step by step. I would target your Integrate & Sum.

I'd be writing all results at once instead of using OpenAppend too.

As a guide, on my Linux 18.04 XUbuntu 12.0 combination (Xeon E5-2690 v4 @ 2.60GHz), the calculation Y1[10., 1.] takes 221.04 seconds.

For Y1[1., 1.] I get lots of underflows. As an example

Timing[Y1[1., #]]& /@ {-1, 0, 1}
During evaluation of In[20]:= General::munfl: Exp[-799.438] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::munfl: Exp[-986.96] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::munfl: Exp[-1194.22] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::stop: Further output of General::munfl will be suppressed during this calculation.
Out[20]= {{226.588,0.247026},{224.858,-0.999189},{224.499,0.247026}}

so the Y1[1., #]& /@ Range[-1., 1, .01] calculation is going to take (at an average of 226 seconds per t) about 12 hours and 40 minutes.

Failing this, I reach for gfortran.