NIntegrate MonteCarlo with phase-space constraints

Question

I wanted to use NIntegrate in Mathematica to perform a multi-dimensional Monte-Carlo integration numerically.

I am a physics student, and I want to calculate probabilty (cross-section) of a high energy scattering event (given their phase-space).

I have simplified the integral to a 3-dim integration to calculate the volume of a sphere (I don't want to do it analytically). The equation is given as:

$$V=\left.\int_{-r}^{+r} ~1~ dx~dy~dz\right|_{x^2+y^2+z^2\leq r^2}$$

So its a homogeneous sphere with a fixed density of 1. The integration phase space is a 3D cube with edge from $-r$ to $+r$. But I am using a constraint to cut certain parts off the phase-space to give it a spherical volume.

The Mathematica code is written as:

(*global*)
r=1;
vol=0.0;

(*integrand*)
f[v1_,v2_,v3_]:=
  Module[{x=v1,y=v2,z=v3},
    dvol=1.0;
    If[x^2+y^2+z^2>=r^2,dvol=0.0];
    dvol
  ]

(*integration routine*)
vol=NIntegrate[
  f[x,y,z],
  {x,-r,+r},
  {y,-r,+r},
  {z,-r,+r},
  Method->{"AdaptiveMonteCarlo","SymbolicProcessing"->0},
  AccuracyGoal->10,
  MaxPoints->10^6,
  MaxRecursion->20
];
vol

which is not what I expected. I coult do a brute force loop which simulates the process:

(*brute loop*)
n=20;
vol=0.0;
For[ix=-n,ix<=+n,ix++,
  For[iy=-n,iy<=+n,iy++,
    For[iz=-n,iz<=+n,iz++,
      x=ix/n; y=iy/n; z=iz/n;
      dvol=(r/n)^3*1.0;
      If[x^2+y^2+z^2<=r^2,vol+=dvol,vol+=0.0]
    ]
  ]
];
vol

and I do get an approximation for the sphere volume. The question is, what should I do to improve on the NIntegrate code above.

Please note that I don't want to change the integrand that much. The integrand should take some input, perform some calculation and spits out a number. The calculation involves multiple conditions (which might not be a simple Heaviside as in the example) and loops, which I enclosed it in a module.

I am guessing this is the parameter tuning for the NIntegrate function. I am not sure what should be the correct setting and I hope someone on SE could help me out. Thanks in advance.

Answer 1

The behaviour you observe is due to symbolic nature of Mathematica. When you write f[x,y,z] in NIntegrate the function f evaluates before x, y, and z receive their numeric values. Since x^2+y^2+z^2>=r^2 (with yet unassigned x, y, z) does not evaluate to True, the dvol=0.0 assignment never happens, and the expression f[x,y,z] evaluates to 1.0.

One way to fix this is to avoid evaluation of f until the numeric arguments are provided:

ClearAll[f];
f[v1_?NumericQ,v2_?NumericQ,v3_?NumericQ]:=
  Module[{x=v1,y=v2,z=v3},
    dvol=1.0;
    If[x^2+y^2+z^2>=r^2,dvol=0.0];
    dvol
  ]

With these modifications, the integral evaluates to the correct answer $4.2$ ($\pm 0.2$). The exact value of the integral is $4\pi/3$, i.e. approximately $4.188790204786390984616857844372670512263$.