0
$\begingroup$

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.

$\endgroup$

1 Answer 1

1
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.