# Is Mathematica really integrating something?

I have a recursive method to get a function, which is already working in other programs, so everything is alright with the method. My problem is when I try to integrate cond[wr] in some limit, like {wr, -1, 1}, it takes so long and I get not from this.

Can you tell me why this NIntegrate is not working?

NIntegrate[cond[wr], {wr, -1, 1}]



Where cond[wr] is obtained from: https://pastebin.com/d5zbYfnH

• I think the problem is simple: Maybe your code is just awfully inefficient. Is see lines like Subscript[b, i] = ... involving Inverse on the right hand side. That somehow tells me, that you might invert the same matrices multiple times. In particular, I found this at several places: Inverse[Energy - H - \[CapitalAlpha] - \[CapitalBeta] - Subscript[CC, i - 1]]. So, I suggest to remove the redundant computations first in order to make your function cond more efficient. – Henrik Schumacher Aug 9 '19 at 21:51
• Moreover, using LinearSolve instead of Inverse should be faster and more accurate. With Inverse, you are likely having numerical instabilities that leads to errors that are more or less randomly distributed. And that might also hinder NIntegrate from converging: NIntergrate supposes that the integrand is more or less smooth; but it has a hard time to deal with a function with basically random function values. – Henrik Schumacher Aug 9 '19 at 21:54
• In which way I could write the coeficients "Subscript[b, i] " in therms of LinearSolve? – Lucas Lopes Aug 10 '19 at 4:41

# Code Refactoring

This is the part that I did not touch, since it is decently efficient and is called only once:

getPandH[wr_] :=
Module[{n = 7, s = 20, η = 10^-3, ϵ = 0, t = 2.75, U, HH, Energy, L},
If[OddQ[n], L = (n - 1)/2;
U = DiagonalMatrix[Flatten[{0, t, 0, Table[{0, 0, 0, t, 0, 0, 0, t}, {i, 1, L}]}], -3];
HH = DiagonalMatrix[Flatten[{ϵ, ϵ, ϵ, ϵ, ϵ, ϵ,Table[{ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ}, {i, 1, L}]}], 0]
+ DiagonalMatrix[Flatten[{t, t, t, t, t, Table[{0, t, t, t, 0, t, t, t}, {i, 1, L}]}], 1]
+ DiagonalMatrix[Flatten[{0, 0, 0, Table[{0, 0, 0, 0, t, 0, 0, 0}, {i, 1, L}]}], 3]
+ DiagonalMatrix[Flatten[{0, 0, Table[{t, 0, 0, 0, 0, 0, 0, 0}, {i, 1, L}]}], 4]
+ DiagonalMatrix[Flatten[{t, Table[{0, 0, 0, 0, 0, 0, 0, t}, {i, 1, L}]}], 5] +
If[n == 1, 0,
DiagonalMatrix[
Flatten[Table[{0, 0, 0, t, 0, 0, 0, 0}, {i, 1, L}]], 6]];
Energy = (wr + I*η)*
DiagonalMatrix[
Flatten[{1, 1, 1, 1, 1, 1,
Table[{1, 1, 1, 1, 1, 1, 1, 1}, {i, 1, L}]}], 0];, L = n/2;
U = DiagonalMatrix[
Flatten[{0, t, 0, 0, 0, 0, t,
Table[{0, 0, 0, t, 0, 0, 0, t}, {i, 1, L - 1}]}], -3];
HH = DiagonalMatrix[
Flatten[{ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ,
Table[{ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ, ϵ}, {i, 1, L - 1}]}],
0] + DiagonalMatrix[
Flatten[{t, t, t, t, t, 0, t, t, t,
Table[{0, t, t, t, 0, t, t, t}, {i, 1, L - 1}]}], 1] +
DiagonalMatrix[
Flatten[{0, 0, t, 0, 0, 0,
Table[{0, 0, 0, 0, 0, 0, 0, 0}, {i, 1, L - 1}]}], 4] +
DiagonalMatrix[
Flatten[{0, 0, 0, 0, 0, 0, 0,
Table[{t, 0, 0, 0, t, 0, 0, 0}, {i, 1, L - 1}]}], 3] +
DiagonalMatrix[
Flatten[{t, 0, 0, 0, 0,
Table[{0, 0, 0, t, 0, 0, 0, t}, {i, 1, L - 1}]}], 5] +
DiagonalMatrix[
Flatten[{0, 0, 0, t,
Table[{0, 0, 0, 0, 0, 0, 0, 0}, {i, 1, L - 1}]}], 6];
Energy = (wr + I*η)*
DiagonalMatrix[
Flatten[{1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
Table[{1, 1, 1, 1, 1, 1, 1, 1}, {i, 1, L - 1}]}], 0];
];
{Energy - (ConjugateTranspose[HH] + HH), U}
];


This is the main code for which I had to invest quite alot of time to refactor because of the extensive use of Subscript which completely obfuscated what was happening there. Also, the algorithm with Subscript required about 10 times more memory than actually needed. In general, Subscript is evil; it does not behave as you expect, and it should be generally avoided. A better option is using downvalues, e.g. a[i] instead of Subscript[a,i]. But even better is not storing data that you don't need.

After refactoring, it occured to me that this is actually a fixed point problem in disguise (hey, you could have told us, right?). So I recast the Table into a call to FixedPoint.

Notice also that variable scoping may be vital for letting this run in parallel.

Clear[cond2];
cond2[wr_, OptionsPattern[{
"Tolerance" -> 1*^-12,
"MaxIterations" -> 20
}]] := Module[{
a, b, c, d, A, B, AA, BB, CC, DD, U, UH, P, L, Energy, α, β, Pαβ, PAB, G0, G2, G11, ΣL, ΣR, ΓL, ΓR, F, residual, TOL, iter, maxiter, X, Y
},
TOL = OptionValue["Tolerance"];
maxiter = OptionValue["MaxIterations"];
{P, U} = getPandH[wr];
UH = ConjugateTranspose[U];

With[{S = LinearSolve[P]},
With[{SU = S[U], SUH = S[UH]},
α = B = UH.SU;
β = A = U.SUH;

Pαβ = P - α - β;
PAB = P - A - B;

a = BB = UH.SUH;
b = AA = U.SU;
c = d = 0.;
CC = DD = 0.;
]
];

F[{a_, b_, c_, d_, AA_, BB_, CC_, DD_}] :=
Block[{Sαβ, S, sa, sb, SAA, SBB},
Sαβ = LinearSolve[Pαβ - c];
S = LinearSolve[PAB - CC];
sa = Sαβ[a];
sb = Sαβ[b];
SAA = S[AA];
SBB = S[BB];
{a.sa, b.sb, c + a.sb + b.sa, d + b.sa, AA.SAA, BB.SBB,
CC + AA.SBB + BB.SAA, DD + BB.SAA}
];

(*loop*)
(*{a,b,c,d,AA,BB,CC,DD}=Nest[F,{a,b,c,d,AA,BB,CC,DD},
s-1];*)
residual = 1.;
iter = 0;
X = {a, b, c, d, AA, BB, CC, DD};
While[residual > TOL && iter < maxiter,
iter++;
Y = F[X];
residual =
Max[Abs[X[[4]] - Y[[4]]]]/Max[Abs[Y[[4]]]] +
Max[Abs[X[[8]] - Y[[8]]]]/Max[Abs[Y[[8]]]];
X = Y;
];
d = X[[4]];
DD = X[[8]];
(*---------------------------------------------------------------------\
*)
G0 = LinearSolve[P - β - d];
G2 = LinearSolve[P - B - DD];

ΣL = U.G0[UH];
ΣR = UH.G2[U];
G11 = LinearSolve[P - ΣL - ΣR];
ΓL = I (ΣL - ConjugateTranspose[ΣL]);
ΓR = I (ΣR - ConjugateTranspose[ΣR]);
Re[Tr[G11[ΓR].G11[ΓL, "J"]]]
]


This is not blazing much faster than the original implementation, though. We get a acceleration facto of about 4, which roughly corresponds to the fact that the original code inverted each matrix four times.

For an argument close to 1.:

cond[0.9] // RepeatedTiming
cond2[0.9] // RepeatedTiming


{0.031, 1.99629}

{0.0069, 1.99629}

For an argument close to 0.:

{0.030, 9.63865*10^-6}

{0.0039, 9.63865*10^-6}

I leave the accuracy checks to the OP.

# Numerical Integration

Let's see what gives NIntegrate such a hard time:

n = 2000;
a = -1;
b = 1;
x = Subdivide[N[a], N[b], n];
vals = ParallelMap[cond2, x];

ListLinePlot[vals, DataRange -> {a, b}]


As we can see, the function has some extremely steep regions and probably thesecause NIntegrate to adaptively refine around these regions forever.

However, we can apply Tai's method to integrate it by hand:

ω = (b - a) ConstantArray[1./n, n + 1];
ω[[1]] *= 0.5;
ω[[-1]] *= 0.5;
ω.vals


2.17057

I tried the same with n = 10000 and the difference of the two results was -1.87783*10^-6, so we may expect that the 5 leading digits are correct.

• First of all, I would like to thanks for your dedication to this subject. I'm learning so much with this post. You got pretty accurate with this reformulation, now the results are calculated with lower computacional costs. How about this Tai's method, what a crazy story, a kind of one rediscovery of Riemann’s sum, very cool. – Lucas Lopes Aug 11 '19 at 4:24
• The original idea were to integrate cond[wr] using Riemann's sum, and plot mu vs sum1/(Te*sum). With kt=kb*Te (kb=8.61*10^-5,Te=300). data[mu_]:=Module[{kt=0.02583,Te=300,sum=0,sum1=0},Do[sum = sum + (f[wr]*0.01); sum1 = sum1 + (f[wr]*(wr - mu)*0.01);, {wr,mu - 50*kt, mu + 50*kt, 0.01}]; sum1/(Te*sum)].  Where f[wr]=E^((wr - mu)/(kt))/(kt*((E^((wr - mu)/(kt)) + 1.)^2))*cond[wr] Then use it as a function itself in a table ListLinePlot[Table[{mu,data[mu]},{mu,-0.4,0.4,0.1}]]. Now with this Tai's method, what would you suggest to implement this? – Lucas Lopes Aug 11 '19 at 4:37
• Using the Riemann's sum above, for n=6 in getPandH[wr_] doing Table[{mu, data[mu]}, {mu, -2.5, 2.5, 0.1}]; I got the results in 684.624 seconds. – Lucas Lopes Aug 11 '19 at 5:04
• Okay, let's look at the body of your function data: There the function f[wr] is called twice for each value of wr. Let's see what f does... ah! It calls cond! Might this be the most expensive computation around? Maybe we should some programming effort to prevent it from beeing called redundantly... – Henrik Schumacher Aug 11 '19 at 12:40
• What about data[mu_] := Module[{kt = 0.02583, Te = 300., sum = 0., sum1 = 0., kb = 8.61*10^-5, a}, Do[ a = E^((wr - mu)/(kt))/(kt ((E^((wr - mu)/(kt)) + 1.)^2)) cond[ wr]; sum = sum + (a 0.01); sum1 = sum1 + (a (wr - mu) 0.01); , {wr, mu - 50*kt, mu + 50*kt, 0.01} ]; sum1/(Te*sum) ]? Then, with n = 6, the computation ParallelMap[data, Subdivide[-2.5, 2.5, 50]] takes about 20 seconds. At least on my Quad Core machine. However, removing the redundance and 4-times acceleration by ParrallelMap can explain only 8-fold speedup. – Henrik Schumacher Aug 11 '19 at 12:43