# How to evaluate a multidimensional NIntegral in a range where certain conditions must be satisfied by the variables

I have a function f which looks similar to this : f[a_, b_, c_] := NIntegrate[sigma[a, b, c], {a, -valueA, valueA}, {b, -valueB, valueB}, {c, -valueC, valueC}]

and wish to integrate it only for valuesa, b, c where another function g fulfills a certain condition : g[a,b,c] >= threshold.

I did try using a boolean in this way: f[a_, b_, c_] := NIntegrate[Boole[g[a, b, c] >= threshold]*sigma[a, b, c], {a, -valueA, valueA}, {b, -valueB, valueB}, {c, -valueC, valueC}]  but I do not get the desired result.

I have also tried to define a Piecewise function for g this way and include it in the integral instead of the Boole: Piecewise[{{g[a,b,c] , g[a,b,c]>= threshold}}]

However, I'm afraid that when using the Piecewise it gets integrated as well, which is not what I wish for. This is just a basic example and in reality I need to pass at least 3 different conditions before I integrate. Looking forward for any tips and help, it's gonna be much appreciated.

tl;dr Trying to numerically integrate a multidmensional integral, and only pass certain values for the variables where conditions a-priori to the integration are fulfilled.

Here's the full integral with prerequisites and values:


(*Transferred energy*)

Tmaxc12[vx_, vy_, vz_, U_, phi_, theta_] :=
0.5*MC12 (vx^2. + vy^2. + vz^2.) + (1 -
Cos[theta])*(Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2] + MC12*vz)*
Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]/MC12 -
Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]*
Sin[theta]*(vx*Cos[phi] + vy*Sin[phi])

(*CONSTANTS DEFINITION*)

Te[U_?NumericQ] := U*e;
\[Beta][U_?NumericQ] := Sqrt[1. - 1./((U/m1) + 1.)^2.];
pe[U_] := Sqrt[Te[U]*(Te[U] + 2.*m*c^2.)/c^2.];
c = 299792458.; (*speed of light*)

m = 9.10938356*10^(-31.);

m1 = 510998.;(*electron mass in eV*)

MC12 = 12.011*1.660539040*10^(-27.);

e = 1.60217662*10^(-19.); (*elementary charge*)
\[HBar] =
1.054571800*10^(-34.); (*reduced Planck constant*)

Zc12 = 6.;

eps = 8.85418*10^(-12. );(*vacuum permittivity*)
(*Velocity \
distributions*)

Pvel[v_?NumericQ, Vfit_?NumericQ] :=
1./Sqrt[2.*Pi*Vfit]*Exp[-v^2./(2.*Vfit)]
(*mean squared velocities for C12*)
VfitxyC12 = 1146080.;
VfitxC12 = VfitxyC12/2.; VfityC12 = VfitxyC12/2.; VfitzC12 = 317000.;
vxvalC12 = Sqrt[VfitxC12]; vyvalC12 = Sqrt[VfityC12]; vzvalC12 =
Sqrt[VfitzC12];

(*cross section*)

k1C12 = ((Zc12 e^2.)/(4. \[Pi] eps 2. m c^2.))^2.;
k2C12 = \[Pi] Zc12 e^2. /(\[HBar] c);
sigmaC12[theta_, U_] :=
k1C12* (1. - \[Beta][U]^2.) /\[Beta][
U]^4.*(Csc[theta/2.])^4.*(1. - \[Beta][U]^2.*Sin[theta/2.]^2. +
k2C12*\[Beta][U]*Sin[theta/2.] (1. - Sin[theta/2.]))*10.^28.;


This is how I defined my region of interest, where Tmax>= 21.14:

region = ImplicitRegion[
Tmaxc12[vx, vy, vz, U, phi, theta]/e >=
21.14, {{vx, -vxvalC12, vxvalC12}, {vy, -vyvalC12,
vyvalC12}, {vz, -vzvalC12, vzvalC12}, {phi, 0, 2 Pi}, {theta, 0,
Pi}}];


and now the integral I was trying to solve :

sigma5D[U_] :=
NIntegrate[
sigmaC12[theta, U]*Sin[theta]*Pvel[vx, VfitxC12]*Pvel[vy, VfityC12]*
Pvel[vz, VfitzC12], {vx, vy, vz, theta, phi} \[Element] region,
sigma5D // Timing


error msg:

The region given at position 1 in DiscretizeRegion[ImplicitRegion[...]] is in dimension 5. DiscretizeRegion only supports dimensions 1 through 3.


after which mathematica crashes and quits the kernel.

• What if you try creating an implicit region reg = ImplicitRegion[g[a, b, c] >= threshold, {a, b, c}] then integrate over that? NIntegrate[sigma[a, b, c], {a, b, c} ∈ reg] Jul 28, 2020 at 11:36
• Hmm.. did not think of it. Will try that. Thanks! Jul 28, 2020 at 12:17
• I tried using the ImplicitRegion method. However i get the following error message, because it is by default trying to discretize the region : DiscretizeRegion[ Implicitregion[...]] is in dimension 5. DiscretizeRegion only supports dimensions 1 through 3 . Any thoughts on how I could bypass this error? The example I referred to is in 3D but I am trying to compute a 5D integral. Jul 28, 2020 at 13:50
• The region of interest is the following :  region = DiscretizeRegion[ ParametricRegion[ Tmaxc12[vx, vy, vz, U, phi, theta] >= threshold, {{vx, -vxval, vxval}, {vy, -vyvalC12, vyvalC12}, {vz, -vzvalC12, vzvalC12}, {phi, 0, 2 Pi}, {theta, 0,Pi}}]]  and the integral I'm trying to solve :  sigma5D[U_] := NIntegrate[ sigmaC12[theta, U]*Sin[theta]*Pvel[vx,VfitxC12]*Pvel[vy,VfityC12]*Pvel[vz, VfitzC12], {vx, vy, vz, theta, phi} [Element] region, Method -> "GlobalAdaptive"]  Jul 28, 2020 at 14:13
• Please put it in the question formatted, and include all values necessary. You have some syntax errors like [Element] and you don't need the DiscretizeRegion. Jul 28, 2020 at 14:16

This is not a complete answer but an extended comment. ImplicitRegion does not like the usage of function Tmaxc12, so we can construct it inline:

region[U_?NumericQ] := ImplicitRegion[
(0.5*MC12 (vx^2. + vy^2. + vz^2.) + (1 -
Cos[theta])*(Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2] + MC12*vz)*
Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]/MC12 -
Sqrt[Te[U]*(Te[U] + 2 m*c^2)/c^2]*
Sin[theta]*(vx*Cos[phi] + vy*Sin[phi]))/e >=
21.14, {{vx, -vxvalC12, vxvalC12}, {vy, -vyvalC12,
vyvalC12}, {vz, -vzvalC12, vzvalC12}, {phi, 0, 2 Pi}, {theta, 0, Pi}}];


Now draw random points from the region, provided U is large enough and the region isn't too 'thin':

pts = RandomPoint[region[100000.], 50000];


Define the integrand:

integrand[U_, {vx_, vy_, vz_, theta_, phi_}] :=
sigmaC12[theta, U]*Sin[theta]*Pvel[vx, VfitxC12]*Pvel[vy, VfityC12]*Pvel[vz, VfitzC12]


We can then look at the values the integrand takes on these points. Notice that they are extremely small in magnitude almost everywhere except at a handful of extreme values.

ListPlot[Sort[integrand[100000., #] & /@ pts], PlotRange -> All]


Re-running the above with different random points will show that negative values and positive values in the tails balance out, while most of the integrand is zero. It's very likely that your integral is zero or so close to zero as to be lost in numerical error.

Trying Monte-Carlo won't settle on any reasonable number for successive runs either:

Mean[integrand[100000.,#]& /@ RandomPoint[region[100000.],50000]]


The following approach will fail too:

With[{reg = region[100000.]},
NIntegrate[
If[RegionMember[reg, {vx, vy, vz, theta, phi}],
integrand[100000., {vx, vy, vz, theta, phi}], 0], {vx, -vxvalC12,
vxvalC12}, {vy, -vyvalC12, vyvalC12}, {vz, -vzvalC12,
vzvalC12}, {phi, 0, 2 Pi}, {theta, 0, Pi}
]]

(* 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 1.2042305021128508322386174756143336864745417085480816121406175838965.954589770191*^645 and 4.3560978955265977448606765353217067111428570569938465058878574724765.954589770191*^643 for the integral and error estimates. *)

(* 1.204230502112851*10^645 *)

• Thank you for your comprehensive answer and analysis to the problem. I figured out why the integrand regions are 0 for most of the part and that is, because the tails of the distributions Pvel x,y,z are not being sampled within higher ranges. It is only for values around the tails of the distributions that TmaxC12 / e  will be higher than 21.14 and here the contributions of  Pvel[vz,VfitzC12]  play the most important role. So the limits should be - 4 vzval, 4vzval and so forth also for vx and vy Jul 29, 2020 at 12:53
• I'm still trying to get around it. The expected calculated integral value should be lying around ~ 0.336. Using the Monte Carlo method is not ideal since for most of the part, the integrand will almost everytime evaluate to 0. GlobalAdaptive seems to be the most well suited method for this purpose, from my understanding so far Jul 29, 2020 at 12:56
• GlobalAdaptive will not work and it's not ideal for higher dimensional integrals. AdaptiveMonteCarlo might be better because it will close in on more interesting areas of higher variance. Jul 29, 2020 at 12:59
• ^ even so, I tried with the wider limits in both the region and the integral, and using AdaptiveMonteCarlo, and running the same experiment in the answer too, but there are still extreme contributions to the integral that lead to unusable results. Jul 29, 2020 at 13:04