# How to speed up a complex function involving integration and matrix inversion?

Here I have a function $$I(V)$$ (current with respect to voltage $$V$$) involving numerical integration of dummy variable $$E$$ for a given $$V$$.

The whole expression of current is given by the following:

$$I(V)=\frac{1}{2\pi}\int_{-\infty}^{\infty} dE \sum_{m=0}^{Nm}\sum_{m'=0}^{Nm'}e^{-E_m/kT}\left[T^{L\rightarrow R}_{m\rightarrow m'}(E,V)f_L(E-E_m;V)(1-f_R(E-E_{m'};V))\\-T^{R\rightarrow L}_{m\rightarrow m'}(E,V)f_L(E-E_m;V)(1-f_R(E-E_{m'};V))\right]$$

where

$$f_{\alpha}(E,V)=(e^{(E-\mu_{\alpha}(V))/kT}+1)^{-1}$$, $$\mu_{\alpha}=(-1)^{\alpha}V/2$$ ($$\alpha=0,1$$ for L and R respectively)

Now comes the bulky part

$$T^{L\rightarrow R}_{m\rightarrow m'}(E,V)=\Gamma^R(E-E_{m'})|G_{m'm}(E,V)|^2\Gamma^R(E-E_{m})$$

with inverse matrix $$G^{-1}_{m'm}(E,V)=F_{m'm}(E,V)=\delta_{m'm}\left[E-\Sigma_L(E-E_{m})-\Sigma_R(E-E_{m})\right]-\sum_{\nu=0}^{N\nu}E_{\nu}C_{m'\nu}C_{m\nu}$$

The other functions like $$\Gamma(E,V),\Sigma(E,V)$$ are usual function defined in the code (see below). The matrix element $$C_{m\nu}$$ is defined as $$C_{m\nu}=e^{-\lambda^2/2\Omega^2}\sum_{i=0}^{\nu}\sum_{j=0}^{m}\delta_{m-j,\nu-i}(-1)^j\left(\frac{\lambda}{\Omega}\right)^{i+j}\frac{1}{i!j!}\sqrt{\frac{m!\nu!}{(m-j)!(v-i)!}}$$.

To evaluate the $$\sum_{m'=0}^{Nm'}T^{L\rightarrow R}_{m\rightarrow m'}(E,V=0)$$, it already takes me a bit long (I plot for $$m=0,1,2$$).

Plot[{\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[0, mf, E, 0]$$\), \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[1, mf, E, 0]$$\), \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[2, mf, E, 0]$$\)}, {E, -0.5, 3}, PlotRange -> All] // AbsoluteTiming


To calculate the current it took me 3 hours to finish. Therefore, I am wondering if there is anything can help speed up this.

Code is below:

Clear["Global"];
\[Lambda] = 0.3;
\[CapitalOmega] = 0.5;
\[Epsilon]0 = 0.5;
\[Beta]1 = 1; \[Beta]2 = 0.2; kT = 0.0259/300*10; (*10K*)
Nm = 7; N\[Nu] = Nm; (*Nm=7 is good enough*)
(*Franck-Condon factors*)

FK = E^(-\[Lambda]^2/(2 \[CapitalOmega]^2)) Table[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 0$$, $$\[Nu]$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 0$$, $$m$$]$$(KroneckerDelta[m - j, \[Nu] - i])$$
\*SuperscriptBox[$$(\(-1$$)\), $$j$$]
\*SuperscriptBox[$$( \*FractionBox[\(\[Lambda]$$, $$\[CapitalOmega]$$])\), $$i + j$$]
\*FractionBox[$$1.$$, $$\(i!$$ $$j!$$\)]
\*SqrtBox[
FractionBox[$$\(m!$$ $$\[Nu]!$$\), $$\(\((m - j)$$!\) $$\((\[Nu] - i)$$!\)\)]]\)\), {m, 0, Nm}, {\[Nu], 0, N\[Nu]}];
(*Eigenenergies of oscillator*)
Em = Table[\[CapitalOmega] (m + 1/2), {m, 0, Nm}];
(*Eigenenergies of shifted oscillator*)
E\[Nu] = Table[\[CapitalOmega] (\[Nu] + 1/
2) + \[Epsilon]0 - \[Lambda]^2/\[CapitalOmega], {\[Nu], 0,
N\[Nu]}];
\[Mu][\[Alpha]_, qV_] := (-1)^\[Alpha] qV/2.;
(*Band width*)
\[CapitalGamma][E_, \[Alpha]_, qV_] :=
If[Abs[(E - \[Mu][\[Alpha], qV])] <=
2 \[Beta]1, \[Beta]2^2/\[Beta]1^2 Sqrt[
4 \[Beta]1^2 - (E - \[Mu][\[Alpha], qV])^2], 0];
(*Self energy*)
\[CapitalSigma][\[Epsilon]_, \[Alpha]_, qV_] := \[Beta]2^2/(
2 \[Beta]1^2)
Piecewise[{{\[Epsilon] - \[Mu][\[Alpha], qV] -
Sqrt[(\[Epsilon] - \[Mu][\[Alpha],
qV])^2 - (2 \[Beta]1)^2], \[Epsilon] > (\[Mu][\[Alpha],
qV] + 2 \[Beta]1)}, {\[Epsilon] - \[Mu][\[Alpha], qV] -
I Sqrt[(2 \[Beta]1)^2 - (\[Epsilon] - \[Mu][\[Alpha], qV])^2],
Abs[\[Epsilon] - \[Mu][\[Alpha], qV]] <=
2 \[Beta]1}, {\[Epsilon] - \[Mu][\[Alpha], qV] +
Sqrt[(\[Epsilon] - \[Mu][\[Alpha],
qV])^2 - (2 \[Beta]1)^2], \[Epsilon] < (\[Mu][\[Alpha],
qV] - 2 \[Beta]1)}}] ;
(*Inverse of dot Green's function*)
F[\[Epsilon]_, V_] :=
Table[KroneckerDelta[mi,
mf] (\[Epsilon] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 0,
V] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 1, V]) - \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$\[Nu] = 0$$, $$N\[Nu]$$]$$E\[Nu][\([$$$$\[Nu] + 1$$$$]$$] FK[$$[$$$$mf + 1, \[Nu] + 1$$$$]$$] FK[$$[$$$$mi + 1, \[Nu] + 1$$$$]$$]\)\), {mi, 0, Nm}, {mf, 0, Nm}];
(*Dot Green's function*)
G[\[Epsilon]_, V_] := Inverse[F[\[Epsilon], V]];
(*Transmission probability from left to right*)
Trl[mi_, mf_, \[Epsilon]_,
V_] := \[CapitalGamma][\[Epsilon] - Em[[mf + 1]], 1, V] Abs[
G[\[Epsilon], V][[mf + 1, mi + 1]]]^2 \[CapitalGamma][\[Epsilon] -
Em[[mi + 1]], 0, V]
(*Transmission probability from right to left *)
Tlr[mi_, mf_, \[Epsilon]_,
V_] := \[CapitalGamma][\[Epsilon] - Em[[mf + 1]], 0, V] Abs[
G[\[Epsilon], V][[mf + 1, mi + 1]]]^2 \[CapitalGamma][\[Epsilon] -
Em[[mi + 1]], 1, V]
(*Fermi distribution *)
f[E_, \[Alpha]_, qV_] := (Exp[(E - \[Mu][\[Alpha], qV])/kT] + 1)^-1;
(* Inelastic scattering tunneling current*)
is[V_] := Re[1/(2 \[Pi]) NIntegrate[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mi = 0$$, $$Nm$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]
\*SuperscriptBox[$$E$$,
FractionBox[$$-Em[\([$$$$mi + 1$$$$]$$]\), $$kT$$]] $$(Trl[mi, mf, \[Epsilon], V] f[\[Epsilon] - Em[\([$$$$mi + 1$$$$]$$], 0,
V] $$(1 - f[\[Epsilon] - Em[\([$$$$mf + 1$$$$]$$], 1, V])\) -
Tlr[mi, mf, \[Epsilon],
V] f[\[Epsilon] - Em[$$[$$$$mi + 1$$$$]$$], 1,
V] $$(1 - f[\[Epsilon] - Em[\([$$$$mf + 1$$$$]$$], 0,
V])\))\)\)\), {\[Epsilon], -\[Infinity], \[Infinity]}]]


Update: I am trying to compute the current involving numerical integration. However, I found when I integrate the function $$T(E,V)$$. It issues some error due to singularity and become slow down. I also want to know if there is any potential improvement can be done on the function is[V].

Here is the improved code

Clear["Global"];
\[Lambda] = 0.3;
\[CapitalOmega] = 0.5;
\[Epsilon]0 = 0.5;
\[Beta]1 = 1; \[Beta]2 = 0.2; kT = 0.0259/300*10; (*10K*)
Nm = 7; N\[Nu] = Nm; (*Nm=7 is good enough*)
(*Franck-Condon factors*)

FK = E^(-\[Lambda]^2/(2 \[CapitalOmega]^2)) Table[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 0$$, $$\[Nu]$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 0$$, $$m$$]$$(KroneckerDelta[m - j, \[Nu] - i])$$
\*SuperscriptBox[$$(\(-1$$)\), $$j$$]
\*SuperscriptBox[$$( \*FractionBox[\(\[Lambda]$$, $$\[CapitalOmega]$$])\), $$i + j$$]
\*FractionBox[$$1.$$, $$\(i!$$ $$j!$$\)]
\*SqrtBox[
FractionBox[$$\(m!$$ $$\[Nu]!$$\), $$\(\((m - j)$$!\) $$\((\[Nu] - i)$$!\)\)]]\)\), {m, 0, Nm}, {\[Nu], 0, N\[Nu]}];
(*Eigenenergies of oscillator*)
Em = Table[\[CapitalOmega] (m + 1/2), {m, 0, Nm}];
(*Eigenenergies of shifted oscillator*)
E\[Nu] = Table[\[CapitalOmega] (\[Nu] + 1/
2) + \[Epsilon]0 - \[Lambda]^2/\[CapitalOmega], {\[Nu], 0,
N\[Nu]}];
\[Mu][\[Alpha]_, qV_] := (-1)^\[Alpha] qV/2.;
(*Band width*)

cExp = Compile[{{x, _Real}}, If[x > -300., Exp[x], 0.],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];
cf = Compile[{{x, _Real}},
If[x < -300., 1., If[x > 300., 0., 1./(1. + Exp[x])]],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];
c\[CapitalGamma]0 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}}, If[
Abs[E - qV/2.] <=
2 \[Beta]1, \[Beta]2^2/\[Beta]1^2 Sqrt[
4 \[Beta]1^2 - (E - qV/2.)^2], 0.], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"];
c\[CapitalGamma]1 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}}, If[
Abs[E + qV/2.] <=
2 \[Beta]1, \[Beta]2^2/\[Beta]1^2 Sqrt[
4 \[Beta]1^2 - (E + qV/2.)^2], 0.], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"];
c\[CapitalSigma]0 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}}, \[Beta]2^2/(2 \[Beta]1^2)
Piecewise[{{E - qV/2 - Sqrt[(E - qV/2)^2 - (2 \[Beta]1)^2],
E > (qV/2 + 2 \[Beta]1)}, {E - qV/2 -
I Sqrt[(2 \[Beta]1)^2 - (E - qV/2)^2],
Abs[E - qV/2] <= 2 \[Beta]1}, {E - qV/2 +
Sqrt[(E - qV/2)^2 - (2 \[Beta]1)^2],
E < (qV/2 - 2 \[Beta]1)}}], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"];
c\[CapitalSigma]1 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}}, \[Beta]2^2/(2 \[Beta]1^2)
Piecewise[{{E + qV/2 - Sqrt[(E + qV/2)^2 - (2 \[Beta]1)^2],
E > (-qV/2 + 2 \[Beta]1)}, {E + qV/2 -
I Sqrt[(2 \[Beta]1)^2 - (E + qV/2)^2],
Abs[E + qV/2] <= 2 \[Beta]1}, {E + qV/2 +
Sqrt[(E + qV/2)^2 - (2 \[Beta]1)^2],
E < (-qV/2 - 2 \[Beta]1)}}], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"];

(*Has to be compute only once.*)
v = cExp[-Em/kT];
A1 = -FK.(E\[Nu] FK\[Transpose]);

ClearAll[Fnew];
Fnew[\[Epsilon]_, V_] :=
Plus[DiagonalMatrix[
Table[\[Epsilon] -
c\[CapitalSigma]0[\[Epsilon] - Em[[mi + 1]],
V, \[Beta]1, \[Beta]2] -
c\[CapitalSigma]1[\[Epsilon] - Em[[mi + 1]],
V, \[Beta]1, \[Beta]2], {mi, 0, Nm}]], A1];
ClearAll[Gnew];
Gnew[\[Epsilon]_, V_] := Inverse[Fnew[\[Epsilon], V]];

(*Transmission probability from left to right*)
Trl[mi_, mf_, \[Epsilon]_, V_] :=
c\[CapitalGamma]1[\[Epsilon] - Em[[mf + 1]],
V, \[Beta]1, \[Beta]2] Abs[
Gnew[\[Epsilon], V][[mf + 1,
mi + 1]]]^2 c\[CapitalGamma]0[\[Epsilon] - Em[[mi + 1]],
V, \[Beta]1, \[Beta]2]
(*Transmission probability from right to left *)
Tlr[mi_, mf_, \[Epsilon]_, V_] :=
c\[CapitalGamma]0[\[Epsilon] - Em[[mf + 1]],
V, \[Beta]1, \[Beta]2] Abs[
Gnew[\[Epsilon], V][[mf + 1,
mi + 1]]]^2 c\[CapitalGamma]1[\[Epsilon] - Em[[mi + 1]],
V, \[Beta]1, \[Beta]2]
(* Inelastic scattering tunneling current*)
is[V_] := Re[1./(2 \[Pi]) NIntegrate[\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mi = 0$$, $$Nm$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]cExp[$$-Em[\([$$$$mi + 1$$$$]$$]\)/
kT] $$(Trl[mi, mf, \[Epsilon], V] cf[ \*FractionBox[\(\[Epsilon] - Em[\([$$$$mi + 1$$$$]$$] -
V/2\), $$kT$$]] $$(1 - cf[ \*FractionBox[\(\[Epsilon] - Em[\([$$$$mf + 1$$$$]$$] +
V/2\), $$kT$$]])\) - Tlr[mi, mf, \[Epsilon], V] cf[
\*FractionBox[$$\[Epsilon] - Em[\([$$$$mi + 1$$$$]$$] +
V/2\), $$kT$$]] $$(1 - cf[ \*FractionBox[\(\[Epsilon] - Em[\([$$$$mf + 1$$$$]$$] -
V/2\), $$kT$$]])\))\)\)\), {\[Epsilon], -\[Infinity], \
\[Infinity]}]]
Plot[{\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[0, mf, E, 0]$$\), \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[1, mf, E, 0]$$\), \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[2, mf, E, 0]$$\)}, {E, -0.5, 3}, PlotRange -> All] // AbsoluteTiming
Plot[{\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]$$Trl[0, mf, E, 0]$$\)}, {E, -0.5, 3}, PlotRange -> All] // AbsoluteTiming
NIntegrate[
Trl[0, 0, \[Epsilon], 0], {\[Epsilon], -2, 2}] // AbsoluteTiming
Plot[is[V], {V, 0, 2.5}] // AbsoluteTiming


Update After implementing above and calculate the function is[V]. It's a bit surprising that there is no improvement on it. See figure below (still took 3 hours):

Phew, this is quite a long and dense code, so I can give only some general hints.

First you should look where the actual bottlenecks are. Please do not expect this from members of this forum.

Anyways, here are some spots to look:

The code throws many General::munfl errors because you feed the exponential functions with many negative values of large magniture. Error handling costs a lot of time. Instead of Exp, you may try to use the function

cExp = Compile[{{x, _Real}},
If[x > -300., Exp[x], 0.],
CompilationTarget -> "C"
];


in order to truncate the results to 0. and to avoid error messages.

You call Trl[mi, mf, \[Epsilon], V] very frequently and each call enforces to reasseble the matrix (call to F) and inverting it (call to G). Don't do that. For given values of \[Epsilon] and V, compute the inverse matrix only once.

Also the function f called way more often than required. You have to call it only once for each pair of {\[Epsilon], mf}, and {\[Epsilon], mi}. So it suffices to call it 4 Nm instead of 4 Nm Nm times...

Moreover, you extensively use KroneckerDelta with double sums and tables. This might not be the biggest problem here, but you should in general avoid that (a general double sum with $$N$$ summands has complexity $$N^2$$ while a double sum with KroneckerDeltas can be rephrase into a summation with only $$N$$ summands).

The resulting funtions seems to be quite smooth. So it might suffice to compute it on a relatively coarse grid and use Interpolate.

Well there are certainly dozens of other suggestions I can make. To give you an idea of order the potential improvement, I am positive that one can get a descent plot of is with just a few seconds, maybe a minute of runtime.

Edit 1 -- remark on function F

A1 = -FK . (E\[Nu] FK\[Transpose]);
Fnew[\[Epsilon]_, V_] := Plus[
DiagonalMatrix[
Table[\[Epsilon] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 0,
V] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 1, V], {mi, 0,
Nm}]
],
A1
]

\[Epsilon] = 1;
V = 1.;
aa = F[\[Epsilon], V]; // RepeatedTiming // First
bb = Fnew[\[Epsilon], V]; // RepeatedTiming // First
Max[Abs[aa - bb]]


0.000559539

0.0000533419

4.44089*10^-16

Further, the diagonal could be generated by a compiled function...

Edit 2 -- remark on G

To highlight that storing intermediate results may work miracles:

ClearAll[Gnew];
Gnew[\[Epsilon]_, V_] := Inverse[Fnew[\[Epsilon], V]];

Do[G[\[Epsilon], V];, {1000}] // AbsoluteTiming // First

fixedG = Gnew[\[Epsilon], V];
Do[fixedG, {1000}]; // AbsoluteTiming // First


0.609069

0.00003

Edit 3 -- On the double sum **

It is actually of the form $$u^T A \,v$$ with a matrix $$A$$ and two vectors $$u$$ and $$v$$. Thus we can exploit fast matrix-vector rountines delivered by Dot.

Some preparations:

cExp = Compile[{{x, _Real}},
If[x > -300., Exp[x], 0.],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
cf = Compile[{{x, _Real}},
If[x < -300.,
1.,
If[x > 300.,
0.,
1./(1. + Exp[x])
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
c\[CapitalGamma]0 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}},
If[Abs[E - qV/2.] <= 2 \[Beta]1,
\[Beta]2^2/\[Beta]1^2 Sqrt[4 \[Beta]1^2 - (E - qV/2.)^2],
0.
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
] ;
c\[CapitalGamma]1 =
Compile[{{E, _Real}, {qV, _Real}, {\[Beta]1, _Real}, {\[Beta]2, \
_Real}},
If[Abs[E + qV/2.] <= 2 \[Beta]1,
\[Beta]2^2/\[Beta]1^2 Sqrt[4 \[Beta]1^2 - (E + qV/2.)^2],
0.
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
] ;

(*Has to be compute only once.*)
v = cExp[-Em/kT];
A1 = -FK . (E\[Nu] FK\[Transpose]);

ClearAll[Fnew];
Fnew[\[Epsilon]_, V_] := Plus[
DiagonalMatrix[
Table[\[Epsilon] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 0, V] - \[CapitalSigma][\[Epsilon] - Em[[mi + 1]], 1, V], {mi, 0, Nm}]
],
A1
];
ClearAll[Gnew];
Gnew[\[Epsilon]_, V_] := Inverse[Fnew[\[Epsilon], V]];


Now compare these:

First@RepeatedTiming[
aa = Quiet[
\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$mi = 0$$, $$Nm$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$mf = 0$$, $$Nm$$]
\*SuperscriptBox[$$E$$,
FractionBox[$$-Em[\([mi + 1]$$]\), $$kT$$]] $$(f[\[Epsilon] - Em[\([mi + 1]$$], 0, V] Trl[mi, mf, \[Epsilon],
V] $$(1 - f[\[Epsilon] - Em[\([mf + 1]$$], 1, V])\) -
f[\[Epsilon] - Em[$$[mi + 1]$$], 1, V] Tlr[mi,
mf, \[Epsilon],
V] $$(1 - f[\[Epsilon] - Em[\([mf + 1]$$], 0, V])\))\)\)\)
];
]

First@RepeatedTiming[
G2 = Abs[Gnew[\[Epsilon], V]]^2;
\[CapitalGamma]0 =
c\[CapitalGamma]0[\[Epsilon] - Em,  V, \[Beta]1, \[Beta]2];
\[CapitalGamma]1 =
c\[CapitalGamma]1[\[Epsilon] - Em,  V, \[Beta]1, \[Beta]2];
f0 = cf[(\[Epsilon] - Em - V/2.)/kT];
f1 = cf[(\[Epsilon] - Em + V/2.)/kT];

bb = Subtract[(v \[CapitalGamma]0 f0).G2.(\[CapitalGamma]1 (1. - f1)), (v f1 \[CapitalGamma]1).G2.(\[CapitalGamma]0 (1. - f0))];
]
(aa - bb)/aa


0.0847156

0.000167522

3.19094*10^-15

This is an over 500-fold speed-up...

• Thanks for your reply. Actually, what I am trying to do is to make the function T(E,V) to be a compiled interpolated function. However, for given E,V this has to be calculated at least once. Jul 29, 2021 at 12:43
• The bottleneck is definitely from the function T(E,V) and matrix inversion. For the first plot, it took me 2 minutes. How to compile such matrix inversion procedure? Jul 29, 2021 at 12:46
• The matrix inversion is already compiled. The point is to not run the inversion more often as needed. Jul 29, 2021 at 12:50
• Thanks a lot. Now it runs at least 10 times faster... Jul 29, 2021 at 13:00
• Plotting the integrand, it looks as if its support were typically quite small. Determining the support and communicating it to NIntegrate should help. Precisionwise, I have a bit of stomach ache because the numbers with which we deal here are really tiny. It might be a good idea to rescale your parameters to convert to "natural units". Jul 30, 2021 at 16:21