# Improve speed of NIntegrate for functions which are sums of complex exponentials

I need to improve the efficiency for numerically evaluating a matrix element, which uses NIntegrate of the several functions (which are all sums of complex exponentials).

Creating the functions from eigenvectors of some matrix Hamiltonian:

Clear[eEigens, eEnergy, eFunc, hEigens, hEnergy, hFunc]

eEigens[β_, γ_, dtoR_, f_] :=
eEigens[β, γ, dtoR, f] = Block[{},
m = 10;
ham = Table[i^2 KroneckerDelta[i, j], {i, -m, m}, {j, -m, m}]
+ Table[(-β*dtoR*(1 - γ) - I f)*
KroneckerDelta[i, j + 1], {i, -m, m}, {j, -m, m}]
+ Table[(-β*dtoR*(1 - γ) + I f)*
KroneckerDelta[i, j - 1], {i, -m, m}, {j, -m, m}]
+ Table[-(1/4) β (1 + γ)*
KroneckerDelta[i, j + 2], {i, -m, m}, {j, -m, m}]
+ Table[-(1/4) β (1 + γ)*
KroneckerDelta[i, j - 2], {i, -m, m}, {j, -m, m}];

eEigs =
SortBy[Transpose[Eigensystem[ham, Method -> "Banded"]], First]
] (* Diagonalize hamltonian and store results *)

eEnergy[β_, γ_, dtoR_, f_, n_] :=
eEnergy[β, γ, dtoR, f, n] =
eEigens[β, γ, dtoR, f][[n + 1]][[1]]

eFunc[β_, γ_, dtoR_, f_, n_] :=
eFunc[β, γ, dtoR, f, n] =
Chop[1./Sqrt[N[2 π]]*
eEigens[β, γ, dtoR, f][[n + 1]][[2]]].Table[
Exp[I*j*φ], {j, -m, m}];

hEigens[μ_, β_, γ_, dtoR_, f_] :=
hEigens[μ, β, γ, dtoR, f] = Block[{},
m = 10;
ham = Table[μ*i^2 KroneckerDelta[i, j], {i, -m, m}, {j, -m,
m}]
+ Table[(β*dtoR*(1 - γ) - I f)*
KroneckerDelta[i, j + 1], {i, -m, m}, {j, -m, m}]
+ Table[(β*dtoR*(1 - γ) + I f)*
KroneckerDelta[i, j - 1], {i, -m, m}, {j, -m, m}]
+ Table[
1/4 β (1 + γ)*KroneckerDelta[i, j + 2], {i, -m,
m}, {j, -m, m}]
+ Table[
1/4 β (1 + γ)*KroneckerDelta[i, j - 2], {i, -m,
m}, {j, -m, m}];
hEigs =
SortBy[Transpose[Eigensystem[ham, Method -> "Banded"]], First]
] (* Diagonalize  hamltonian *)

hEnergy[μ_, β_, γ_, dtoR_, f_, n_] :=
hEnergy[μ, β, γ, dtoR, f, n] =
hEigens[μ, β, γ, dtoR, f][[n + 1]][[
1]] (* Read stored results to give eigenvalues *)

hFunc[μ_, β_, γ_, dtoR_, f_, n_] :=
hFunc[μ, β, γ, dtoR, f, n] =
Chop[1./Sqrt[N[2 π]]*
hEigens[μ, β, γ, dtoR, f][[n + 1]][[
2]]].Table[Exp[I*j*φ], {j, -m, m}];


So eFunc and hFunc are sums of complex exponentials for each eEigs and hEigs respectively, with coefficients given by the corresponding eigenvectors eEigens and hEigens.

I then create the following function, which will then be evaluated 100's of times to populate a matrix.

integral[μ_?NumericQ, β_?NumericQ, γ_?NumericQ,
dtoR_?NumericQ, f_?NumericQ, n0_, M0_, n_, M_] := NIntegrate[
eFunc[β, γ, dtoR, f, n]*
hFunc[μ, β, γ, dtoR, f, M]*
Conjugate[eFunc[β, γ, dtoR, f, n0]]*
Conjugate[hFunc[μ, β, γ, dtoR, f, M0]
], {φ, -π, π}, AccuracyGoal -> 10,
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", "SymbolicProcessing" -> 0}
]


This takes about 0.02 seconds on my laptop:

integral[0.4, 2., 1.02, 20., 0., 1, 1, 2, 2] // AbsoluteTiming
(* {0.020465, -0.00181996} *)


Is there a quicker way of doing this?

I then create a final matrix for example:

finalmatrix[μ_, β_, γ_, dtoR_, f_, v0_, x_, max_] :=
Rationalize[
Block[{},
g = Table[(-N[2 π]*v0)/(
eEnergy[β, γ, dtoR, f, jj] +
hEnergy[μ, β, γ, dtoR, f, kk] - x)*
integral[μ, β, γ, dtoR, f, ii, ff, jj, kk],
{ii, 0, max}, {ff, 0, max}, {jj, 0, max}, {kk, 0, max}];
h = ArrayReshape[g, {(max + 1)^2, (max + 1)^2}] -
IdentityMatrix[(max + 1)^2]
]
, 0];


which takes far too long for the purposes of subsequent calculations

finalmatrix[0.4, 2., 1.02, 20., 0., -0.2, x, 9]; // AbsoluteTiming
(* {103.234, Null} *)


Any help will be hugely appreciated!

• You could cut the time by running ParallelTable assuming you have a multi-core machine. Say off-load {kk,0,max/2} to one core and {kk,max/2+1,max} to the other. – Dominic Aug 7 '19 at 21:11
• If you can reformulate your computations as an integral over a list/matrix of functions you can apply the approach described here. That might give you a fairly significant speed-up. – Anton Antonov Aug 9 '19 at 13:00