# NIntegral takes long time

I am doing an Nintegral on Mathmetica. F[T,px,py,pz] is integrand with long expression. It has poles but with negative and positive infinity from two sides. Evaluate it at for example Eq[0.75], Global adaptive give this message and huge error

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one
of the following: singularity, value of the integration is 0, highly
oscillatory integrand, or WorkingPrecision too small.


If I use local adaptive, the integral takes long time to evaluate only even if I s set PrecisionGoal -> 2. How can I do this integral properly in MMA?

d = 0.25; U1 = 0.52; U2 = 0.3;
UX[px_, py_, pz_] = 2 U1/3 (Cos[px] + Cos[py] + Cos[pz]);
UY[px_, py_, pz_] = 2 U2/3 (Cos[px] + Cos[py] - Cos[pz]);
F[T_, px_, py_, pz_] =
With[{Ux = UX[px, py, pz],
Uy = UY[px, py, pz]}, -(1/(
4 d T)) ((2 E^((2 d)/T) T)/(-1 + E^((2 d)/T)) + (
2 T (Sqrt[d^2] Ux + d Uy - Ux Uy Tanh[d/T]))/((-1 + E^((2 d)/
T)) (d (Ux + Uy) -
Ux Uy Tanh[d/T])) + (Sech[d/
T]^2 (d^3 (Ux + Uy)^2 - 3 d^2 Ux Uy (Ux + Uy) Tanh[d/T] +
2 Sqrt[d^2] Ux^2 Uy^2 Tanh[d/T]^2) +
T (-2 ((d^2)^(3/2) Ux + d^3 Uy) -
d^2 (Ux^2 - 2 Ux Uy - Uy^2) Tanh[d/T] +
d Ux Uy (3 Ux + Uy) Tanh[d/T]^2 -
2 Ux^2 Uy^2 Tanh[d/T]^3))/((-1 + E^((
2 Sqrt[d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(
Sqrt[d^2] T))) Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])] (d (Ux + Uy) -
Ux Uy Tanh[d/T])) + (E^((
2 Sqrt[d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(
Sqrt[d^2]
T)) (Sech[d/
T]^2 (d^3 (Ux + Uy)^2 -
3 d^2 Ux Uy (Ux + Uy) Tanh[d/T] +
2 Sqrt[d^2] Ux^2 Uy^2 Tanh[d/T]^2) +
T (-2 ((d^2)^(3/2) Ux + d^3 Uy) -
d^2 (Ux^2 - 2 Ux Uy - Uy^2) Tanh[d/T] +
d Ux Uy (3 Ux + Uy) Tanh[d/T]^2 -
2 Ux^2 Uy^2 Tanh[d/T]^3)))/((-1 + E^((
2 Sqrt[d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(
Sqrt[d^2] T))) Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])] (d (Ux + Uy) -
Ux Uy Tanh[d/T])))];
Eq2[T_] :=
NIntegrate[
F[T, px, py,
pz], {px, -\[Pi], \[Pi]}, {py, -\[Pi], \[Pi]}, {pz, -\[Pi], \
\[Pi]}];

• Simplify can make your integrand 2/3 the size in a handful of seconds. Sometimes that can speed up NIntegrate, sometimes not. BUT every time I try to evaluate your code it bails telling me that it was non-numeric. Can you scrape-n-paste exactly that post back into a new notebook, evaluate just that and see if it works for you? And if not then figure out what happened and try to fix that?
– Bill
May 11, 2018 at 20:46
• @Bill corrected it
– p.s
May 11, 2018 at 21:07
• adding Method->"AdaptiveMonteCarlo" to Nintegrate and then Eq2[1] gives an answer of just about 1000 after about fifteen seconds, even without using Simplify to speed this up, but with all the poles, is this really an accurate answer or not?
– Bill
May 11, 2018 at 21:52
• Yes, I need more accurate answer and it will be better If I could compare Montecarlo with other method.
– p.s
May 11, 2018 at 22:08
• More accurate than what? To compare methods you could try each one and see what results you get. You should know that random processes, like MonteCarlo and AdaptiveMonteCarlo will have different results every time they run because they randomly choose different points every time. I urge you to do the comparisons and then edit your original post to show how they compare. And you might see how much Simplify changes the results. You might see if giving Simplify information about the domain of each variable makes a difference. So you have experiments to try. FullSimplify is slow but might help
– Bill
May 11, 2018 at 23:41

Analyse the integrand F to learn a lot, how to do integration.

d = 25/100; U1 = 52/100; U2 = 3/10;
UX[px_, py_, pz_] = 2 U1/3 (Cos[px] + Cos[py] + Cos[pz])
UY[px_, py_, pz_] = 2 U2/3 (Cos[px] + Cos[py] - Cos[pz]);

F[T_, px_, py_, pz_] =
Simplify[With[{Ux = UX[px, py, pz],
Uy = UY[px, py,
pz]}, -(1/(4 d T)) ((2 E^((2 d)/T) T)/(-1 +
E^((2 d)/T)) + (2 T (Sqrt[d^2] Ux + d Uy -
Ux Uy Tanh[d/T]))/((-1 + E^((2 d)/T)) (d (Ux + Uy) -
Ux Uy Tanh[d/T])) + (Sech[d/T]^2 (d^3 (Ux + Uy)^2 -
3 d^2 Ux Uy (Ux + Uy) Tanh[d/T] +
2 Sqrt[d^2] Ux^2 Uy^2 Tanh[d/T]^2) +
T (-2 ((d^2)^(3/2) Ux + d^3 Uy) -
d^2 (Ux^2 - 2 Ux Uy - Uy^2) Tanh[d/T] +
d Ux Uy (3 Ux + Uy) Tanh[d/T]^2 -
2 Ux^2 Uy^2 Tanh[d/T]^3))/((-1 +
E^((2 Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(Sqrt[
d^2] T))) Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])] (d (Ux + Uy) -
Ux Uy Tanh[
d/T])) + (E^((2 Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(Sqrt[
d^2] T)) (Sech[d/T]^2 (d^3 (Ux + Uy)^2 -
3 d^2 Ux Uy (Ux + Uy) Tanh[d/T] +
2 Sqrt[d^2] Ux^2 Uy^2 Tanh[d/T]^2) +
T (-2 ((d^2)^(3/2) Ux + d^3 Uy) -
d^2 (Ux^2 - 2 Ux Uy - Uy^2) Tanh[d/T] +
d Ux Uy (3 Ux + Uy) Tanh[d/T]^2 -
2 Ux^2 Uy^2 Tanh[d/T]^3)))/((-1 +
E^((2 Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])])/(Sqrt[
d^2] T))) Sqrt[
d^2 (d - Ux Tanh[d/T]) (d - Uy Tanh[d/T])] (d (Ux + Uy) -
Ux Uy Tanh[d/T])))], -Pi < px < Pi && -Pi < py <
Pi && -Pi < pz < Pi];


Maximize the functions for some T to see Maximum is at {px,py,pz}=={0,0,0}.

NMaximize[F[3, px, py, pz], {px, py, pz}]

(*   {12.7193, {px -> 7.88861*10^-31, py -> 0., pz -> 0.}}   *)


Plot shows, there are a lot of singularities. Look, where 1/F ==0

Plot[F[T, 0, 0, 0], {T, 0, 2}];

LogLinearPlot[{0, 1/F[T, 0, 0, 0]}, {T, 10^-2, 10}]


The first singularity is at

Solve[1/F[T, 0, 0, 0] == 0 && T > 1/2, T, Reals]

(*   {{T -> 1/(4 ArcTanh[25/104])}} // N  == {{T -> 1.01965}}   *)


Looking for Maximum makes no problems for T > 1.01965, but for T below you get problems because of the singularities.

Plot[NMaximize[F[T, px, py, pz], {px, py, pz}][[1]], {T, 1.02, 10}];

NMaximize[F[1, px, py, pz], {px, py, pz}]

(*   NMaximize::nrnum: The function value 182711+2.93397*10^-9 I
is not a real number at {px,py,pz} = {0.195664,-0.264406,0.0560865}.    >>  *)


Edit 3

See the hypersurface where singularities occur (here for T==3/4, excluding artificial effects at F==0)

ContourPlot3D[
Evaluate[1/F[3/4, px, py, pz]], {px, -Pi, Pi}, {py, 0, Pi}, {pz,  -Pi,
Pi}, Contours -> {0},
RegionFunction ->
Function[{px, py, pz},
1 < Abs[F[3/4, px, py, pz]] || Abs[F[3/4, px, py, pz]] < -1]]


That means, for T > 1/(4 ArcTanh[25/104]), you easily get the integral, but for T below, you have to take the singularities into acount. You have to use Method->"PrincipalValue" in NIntegrate and explicitly give the positions of these singulatities, which is quite complex in this case. (Have a look at the documentation for NIntegrate)

Eq2[T_?NumericQ] := NIntegrate[
F[T, px, py,
pz], {px, -\[Pi], \[Pi]}, {py, -\[Pi], \[Pi]}, {pz, -\[Pi], \
\[Pi]}];

tab1 = Table[Eq2[1.02*i], {i, 1, 60}];

ListLogPlot[tab1, PlotRange -> All]


Edit 4

To get the singularity, look, where numerator of 1/F is zero

ff = Numerator[
1/F[3/4, px, py, pz] // Together //
Simplify[#, 0 < px < Pi && 0 < py < Pi && 0 < pz < Pi] &]


This implies, the expression under the Sqrt to be zero.

Solve[(375 - 20 (41 Cos[px] + 41 Cos[py] + 11 Cos[pz]) Tanh[1/3] +
416 (Cos[px]^2 + 2 Cos[px] Cos[py] + Cos[py]^2 - Cos[pz]^2) Tanh[
1/3]^2) == 0 && 0 < px < Pi && 0 < py < Pi && 0 < pz < Pi, pz]

(*   {{pz -> ConditionalExpression[
ArcCos[1/104 (-104 Cos[px] - 104 Cos[py] + 75 Coth[1/3])],
0 < px < ArcCos[1/104 (-208 + 75 Coth[1/3])] &&
0 < py < ArcCos[1/104 (-104 - 104 Cos[px] + 75 Coth[1/3])]]}}   *)


This defines the z position of the singularity in an restricted area of px and py.

For T == .75, I found singularity positiion beeing an ellipsoid of px,py,pz. For given px,py, singularity is at pz == ArcCos[-Cos[px] - Cos[py] + 75/104 Coth[1/3]].

Since the integrand is totally symmetric

F[3/4, px, py, pz] == F[3/4, -px, -py, -pz] // Simplify

(*   True   *)


Error correction

Since Method -> "PrincipalValue" can only work in areas with a singularity (otherwise ArcCos[-Cos[px] - Cos[py] + 75/104 Coth[1/3]] is imginary, you have to do different integration for the areas with and without singularity.

integrate from 0 to Pi and multiply by factor 8. Method -> "PrincipalValue" works only for 1 dimensional integrals, therefore split integration

(int1 = 8*NIntegrate[
NIntegrate[
F[3/4, px, py, pz], {pz, 0,
ArcCos[-Cos[px] - Cos[py] + 75/104 Coth[1/3]], Pi},
Method -> "PrincipalValue"], {px, 0,
ArcCos[1/104 (-208 + 75 Coth[1/3])]}, {py, 0,
ArcCos[1/104 (-104 - 104 Cos[px] + 75 Coth[1/3])]}]) // Timing

(*   {402.844, -202.995 - 1.11989*10^-12 I}   *)

(int2 = 8*NIntegrate[
F[3/4, px, py, pz], {pz, 0, Pi}, {px, 0,
ArcCos[1/104 (-208 + 75 Coth[1/3])]}, {py,
ArcCos[1/104 (-104 - 104 Cos[px] + 75 Coth[1/3])],
Pi}]) // Timing

(*   {1.312, 499.792}   *)

(int3 = 8*NIntegrate[
F[3/4, px, py, pz], {pz, 0, Pi}, {px,
ArcCos[1/104 (-208 + 75 Coth[1/3])], Pi}, {py, 0,
Pi}]) // Timing

(*   {5.5, 58.6319}   *)

int = int1 + int2 + int3

(*   355.429 - 1.11989*10^-12 I   *)


To get the general z position of the singularity dependend on T needs some further effort.

• For me great answer (+1) May 12, 2018 at 9:53
• How do you find the singularity point pz == ArcCos[-Cos[px] - Cos[py] + 75/104 Coth[1/3]]?
– p.s
May 12, 2018 at 14:43
• Look where ff = Numerator[ 1/F[3/4, px, py, pz] // Together // Simplify[#, 0 < px < Pi && 0 < py < Pi && 0 < pz < Pi] &]  is zero. You have two terms Sqrt[...]* (-1+Exp[2 Sqrt[...]] . with the same Sqrt-term. Solve[Sqrt[ ...] == 0 && 0 < px < Pi && 0 < py < Pi && 0 < pz < Pi, pz]  to get pz -> ArcCos[-Cos[px] - Cos[py] + 75/104 Coth[1/3]]  together with conditions for px and py. May 12, 2018 at 20:10
• Thanks for the great answer. Especially on the part of singularity point. For T<1/4/ ArcCoth[12/5], I find there is one more singular point. I am running int1 and get error message" The integrand has evaluated to Overflow, Indeterminate, or Infinity for all \ sampling points in the region with boundaries ". Is it normal? Also before you split the integral in to three region, the integral is pretty fast, Why now it take much longer time?
– p.s
May 14, 2018 at 16:16
• For lower T, you get more and more singularities, which you all have to take into acount in NIntegrate. See for T ==1/2 ContourPlot3D[ Evaluate[1/F[1/2, px, py, pz]], {px, -Pi, Pi}, {py, 0, Pi}, {pz, -Pi, Pi}, Contours -> {0}, RegionFunction -> Function[{px, py, pz}, 2 < Abs[F[1/2, px, py, pz]] || Abs[F[1/2, px, py, pz]] < -2]]  Looking for the singularities with ff = Numerator[ 1/F[1/2, px, py, pz] // Together // Simplify[#, 0 < px < Pi && 0 < py < Pi && 0 < pz < Pi] &]  yields now two solutions Do Solve. Integration without singulaities is much faster. May 14, 2018 at 18:20