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

Question

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}, 
   Method -> {"GlobalAdaptive"}]]

(*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}}

``` 

Answer 1

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

Answer 2

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} *)

Answer 3

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   *)