# Convergence and efficiency of 5-dimensional numerical integral using “GlobalAdaptive”

We have a 5-dimensional integral arising as a solution to a physics problem, where two dimensions of the integration space cover the spherical and azimuthal angle, and three dimensions are needed to account for Gaussians in three Cartesian directions.

We are using Mathematica's NIntegrate with "GlobalAdaptive", which seems the only way to get any sensible values. The Boolean constraint below is mandatory (we have also tried Heaviside), and we have stripped the physical constants from the expressions below to focus on the problem. Numerically speaking, the problem is that the integrand has a complicated functional dependence (determined by the physics), and the volume of the integration space is very large.

With the settings below, by stepwise increasing the integration limits for the Gaussians, we are able to (probably) reach convergence, but the calculation throws warnings and completely fails for i = 5 (corresponding to five standard deviations from the mean) – before crashing, Mathematica uses 75 GB of memory on a very powerful desktop workstation.

Are there any ways to improve convergence and computational efficiency of the integral below? We need reliable integral values for MANY input parameters, and in the real physical situation the integration is even significantly slower than in this simplified example.

Many thanks!

beta[u_] := Sqrt[1 - 1/(u + 1)^2];

(*integrand*)
sigma[theta_, u_] :=  (1 - beta[u]^2) /
beta[u]^4*(Csc[theta/2])^4*(1 - beta[u]^2*Sin[theta/2]^2 +
beta[u]*Sin[theta/2] (1 - Sin[theta/2]));

(*gaussian distributions*)
gauss[v_, sig_] := 1/(sig*Sqrt[2*Pi])*Exp[-v^2/(2*sig)];

(*function for the boole*)
Tmax[vx_, vy_, vz_, U_, phi_,
theta_] := (1/2)*(vx^2 + vy^2 + vz^2) + (1 -
Cos[theta])*(Sqrt[U*(U + 2)] + vz)*Sqrt[U*(U + 2)] -
Sqrt[U*(U + 2)]*Sin[theta]*(vx*Cos[phi] + vy*Sin[phi])

(*integral*)
fivedim[u_, i_] :=
AbsoluteTiming[
NIntegrate[
Boole[Tmax[cx, cy, cz, u, phi, theta] >= 15]*sigma[theta, u]*
Sin[theta]*gauss[cx, 1]*gauss[cy, 1]*gauss[cz, 1], {theta, 0,
Pi}, {phi, 0, 2 Pi}, {cx, -i, i}, {cy, -i, i}, {cz, -i, i},

(*calculate integral to estimate convergence limit*)
limit = {1, 2, 3, 4};
values = {};
For[i = 1, i <= Length[limit], i++,
values = Append[values, fivedim[100000, i]]]
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.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.5581961194518056 and 0.009589824655407775 for the integral and error estimates.

General::stop: Further output of NIntegrate::eincr will be suppressed during this calculation.


We can see from the results that the value of the integral does not fully converge (see {time [s], integral value} values below).

In[97]:= values

Out[97]= {{8.81786, 0.558196}, {32.2968, 1.64615}, {7.84894,
2.00055}, {10.6271, 2.02892}}

$$$$

• The Append is inefficient, you can just use Array or Table, also you repeatedly calculate the same expressions, you could use With to speed this up. Even so it still has the same issues. See my suggestions here pastebin.com/xFuTpp9f esp. With, but I don't think there's much you can do. Jun 4, 2021 at 14:21
• I agree with flinty that I doubt there are simple tricks that will fix this for you. Have you read this? Maybe a Levin rule is suitable given the periodic functions in your integrand. Jun 4, 2021 at 14:32
• Changing to fivedim[u_, i_] := AbsoluteTiming[ NIntegrate[ Boole[Tmax[cx, cy, cz, u, phi, theta] >= 15]*sigma[theta, u]* Sin[theta]*gauss[cx, 1]*gauss[cy, 1]*gauss[cz, 1], {theta, 0, Pi}, {phi, 0, 2 Pi}, {cx, -i, i}, {cy, -i, i}, {cz, -i, i}, Method -> {"LocalAdaptive"}, AccuracyGoal -> 6, PrecisionGoal -> 6]]; Table[fivedim[100000, i], {i, 1, 5}], I obtain {{5.06368, 1.60554*10^-6}, {5.0857, 0.0000155633}, {6.33888, 0.000208747}, {7.89441, 0.000525671}, {5.09002, 1.58715*10^-6}}.  Jun 4, 2021 at 14:54
• Changing to fivedim[u_, i_] := AbsoluteTiming[ NIntegrate[ Boole[Tmax[cx, cy, cz, u, phi, theta] >= 15]*sigma[theta, u]* Sin[theta]*gauss[cx, 1]*gauss[cy, 1]*gauss[cz, 1], {theta, 0, Pi}, {phi, 0, 2 Pi}, {cx, -i, i}, {cy, -i, i}, {cz, -i, i}, Method -> {"MonteCarlo"}, AccuracyGoal -> 3, PrecisionGoal -> 3]];Table[fivedim[100000, i], {i, 1, 5}], I obtain {{0.0068733, 5.8473*10^-6}, {0.008582, 0.0000275885}, {0.0063151, 0.0000463867}, {0.0068139, 6.62425*10^-7}, {0.0099911, 1.76042*10^-7}}. I think Method->"GlobalAdaptive" is not a proper method here. Jun 4, 2021 at 15:21
• I tried a bit to work with your integral. One interesting thing I found: if you use Expand[] over your whole expression to be integrated, then you end up with 20 sumed-terms. If you integrate with i<5 you'll see, that most of them are insignificant, and there are often pairwise effectively canceling terms. Only one single term does actually contribute to the first 3-4 significant digits. So one could concentrate on just using this single term. However, thats an empirical result. Jun 4, 2021 at 15:22

Regarding the integration speed: changing the method from

Method-> {"GlobalAdaptive"} to Method {"GlobalAdaptive","SymbolicProcessing->0"}


did the trick. The time spent integrating is more than halved without losing any accuracy.

{{3.88798, 0.558196}, {2.88294, 1.69157}, {2.82454,
2.00055}, {3.13186, 2.03672}}


Will continue to look into speeding this further and will update if I find an even better solution.

EDIT: adding "." (e.g 1./2. instead of just 1/2) after the constants in the definitions further improved the speed by about 10% to 15% without affecting the results.

Regarding convergence: still haven't found a way to make the integral converge for increasing limits

Below are some suggestions to explore the computations described in the question.

Rewriting the integration function to take options:

Clear[fivedim2];
fivedim2[u_, i_, opts : OptionsPattern[]] :=
NIntegrate[
Boole[Tmax[cx, cy, cz, u, phi, theta] >= 15]*sigma[theta, u]*
Sin[theta]*gauss[cx, 1]*gauss[cy, 1]*gauss[cz, 1], {theta, 0,
Pi}, {phi, 0, 2 Pi}, {cx, -i, i}, {cy, -i, i}, {cz, -i, i}, opts];


Calls with different settings for precision goal, accuracy, and method:

lsRes =
Map[(
res =
AbsoluteTiming@
fivedim2[
100000, #,
PrecisionGoal -> 6, AccuracyGoal -> 6,
MinRecursion -> 3,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None,
"SymbolicProcessing" -> 0, "MaxErrorIncreases" -> 10000}
];
Print[i, " ", res];
res
) &,
Range[6]
]

(* {{17.5683, 0.130449}, {16.0355, 0.405454},
{17.8186, 0.509822}, {14.0988, 0.493829},
{14.2966, 0.651644}, {15.0894, 0.816988}} *)


Remark: Note that above we have results for i being 5 and 6.

lsRes2 =
Map[(
res =
AbsoluteTiming@
fivedim2[
100000, #,
PrecisionGoal -> 6, AccuracyGoal -> 6,
MinRecursion -> 2,
Method -> {"GlobalAdaptive", "SingularityHandler" -> None,
"SymbolicProcessing" -> 0, "MaxErrorIncreases" -> 10000,
Method -> {"ClenshawCurtisRule", "Points" -> 3}}
];
Print[i, " ", res];
res
) &,
Range[4]
]

(* {35.571, 4.85967*10^-7}, {873.671, 0.0700608}, {853.95, 0.086318}, {985.935,  0.123664} *)

• thank you for your suggestion! This is very helpful for trying different options. However, note that the values of your integral are not converging – i controls the size of the integration volume, so we should be reaching a fixed value as it increases. I was further able to calculate lsRes for i = 7, which gave 0.788444, so while somewhat close to the i = 6 value, still not converged. Higher i are not possible to use on my system. There is a further question of why your integral gives different values than the original one..? Jun 22, 2021 at 9:38

Boole expression for Tmax shows to be strongly dependent on theta. Finding minimal value for Tmax depending on theta (thetamin) and integrating from it, lets NIntegrate work fine without problems.

(Here i regard only u=100000 and i =1)

beta[u_] =
Sqrt[1 - 1/(u + 1)^2] // Together //
PowerExpand[#, Assumptions -> u > 0] &;
sigma[theta_,
u_] = (1 - beta[u]^2)/
beta[u]^4*(Csc[theta/2])^4*(1 - beta[u]^2*Sin[theta/2]^2 +
beta[u]*Sin[theta/2] (1 - Sin[theta/2])) // Together //
FullSimplify[#, Assumptions -> {u > 0, 0 < theta < Pi}] &;
gauss[v_, sig_] = 1/(sig*Sqrt[2*Pi])*Exp[-v^2/(2*sig)];
Tmax[vx_, vy_, vz_, U_, phi_,
theta_] = (1/2)*(vx^2 + vy^2 + vz^2) + (1 -
Cos[theta])*(Sqrt[U*(U + 2)] + vz)*Sqrt[U*(U + 2)] -
Sqrt[U*(U + 2)]*Sin[theta]*(vx*Cos[phi] + vy*Sin[phi]) // Expand

mi[theta_] :=
MinValue[{Tmax[cx, cy, cz, 100000, phi,
theta], -1 <= cx <= 1 && -1 <= cy <= 1 && -1 <= cz <= 1 &&
0 <= phi <= 2 Pi}, {cx, cy, cz, phi}]

LogLogPlot[{15, mi[theta]}, {theta, 10^-10, Pi}, PlotPoints -> 21,
MaxRecursion -> 2]


(Since MinValue generaly knows the rules for given theta values, but checking routines are not mighty enough to check it for general theta, do the trick to give an exotic numerical value in the expected range and reinsert theta.)

mmi[theta_] = mi[EulerGamma/10000] /. EulerGamma -> 10000 theta//Simplify

(*   5000100001 - 5000100000 Cos[theta]^2 +
200 Sqrt[250005]
Sin[theta] (Cos[2 ArcTan[1 + Sqrt[2]]] -
Sin[2 ArcTan[1 + Sqrt[2]]])   *)

{thetamin =
theta /. First@
Solve[mmi[theta] == 15 && 5/100000 < theta < 10/100000, theta],
thetamin // N}

(*   {2 ArcTan[
Root[7 - 22857599972 #1^2 + 14286285674284800042 #1^4 -
22857599972 #1^6 + 7 #1^8 &, 6]], 0.0000689137}   *)

Plot[{15, mi[theta]}, {theta, 0, 1/10000}, PlotPoints -> 21,
MaxRecursion -> 2, GridLines -> {{thetamin}, Automatic}]


NIntegrate[
Boole[Tmax[cx, cy, cz, 100000, phi, theta] >= 15]*
sigma[theta, 100000]*Sin[theta]*gauss[cx, 1]*gauss[cy, 1]*
gauss[cz, 1], {theta, thetamin, Pi}, {phi, 0, 2 Pi}, {cx, -1,
1}, {cy, -1, 1}, {cz, -1, 1}]

(*   0.336783   *)