# NIntegrate MonteCarlo with phase-space constraints

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},
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.

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$$.

• Thank you very much. This trick is very useful for many numerical integration problems. Mar 7, 2023 at 6:23