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

Question

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, 
  Method -> "GlobalAdaptive"]
sigma5D[100000] // 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.

Answer 1

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.20423050211285083223861747561433368647454170854808161214061758389`65.954589770191*^645 and 4.35609789552659774486067653532170671114285705699384650588785747247`65.954589770191*^643 for the integral and error estimates. *)

(* 1.204230502112851*10^645 *)