Numerical 3D Integration in Mathematica [closed]

My question in short is: Is it possible to make numerical 3D integration in mathematica efficiently? Or am I gonna have to move this to C++? I feel this problem is way too trivial to go to C++, so please help.

I tried to find this in Mathematica with built-in functions and I couldn't. So I built my own 3D integrator, but the accuracy isn't so good unless I subdivide my volume into a billion pieces, which will be computationally very expensive.

The problem I'm trying to solve is the following:

I wanna integrate over a 3D object that is defined mathematically. Say it's a sphere. The 3D analytical integral will be:

$$\int_{-R+x_{0}}^{+R+x_{0}}\int_{y_{0}-\sqrt{R^{2}-\left(x-x_{0}\right)^{2}}}^{y_{0}+\sqrt{R^{2}-\left(x-x_{0}\right)^{2}}}\int_{z_{0}-\sqrt{R^{2}-\left(x-x_{0}\right)^{2}-\left(y-y_{0}\right)^{2}}}^{z_{0}+\sqrt{R^{2}-\left(x-x_{0}\right)^{2}-\left(y-y_{0}\right)^{2}}}f(x,y,z)dxdydz$$

where $(x_0, y_0, z_0)$ is the center point of the sphere; and $R$ is its radius. How can I test my integrator's accuracy? I can simply set $f(x,y,z)=1$, and then the integration will result in the volume of the sphere. Comparison with that is very easy with $V=\frac{4}{3}\pi R^3$.

With my numerical integrator, I was able to achieve 0.4% accuracy by comparing with the volume of a sphere. I'm not happy about that, and I need a better solution (and btw, it's very slow even after compiling and parallelizing it with the subdivision volume provided in the code below).

What parameters do I expect to give to the 3D numerical integrator?

I expect to give it the function $f(x,y,z)$ and the boundaries of the integral, which will define my integration volume. Then the integrator should do the job.

My numerical integrator is the following. The algorithm is that the volume is first subdivided into cuboids with dimensions {SubLengthX,SubLengthY,SubLengthZ}, and then the function is evaluated for this cuboid and summed after compiling it.

Notice the parts y+Mod[y2-y1,SubLengthY]/2, which are supposed to shift the mesh to be centered in that axis using the remainder. This also wasn't so pleasant as the plotting the mesh didn't look 100% symmetric for a sphere. That's a separate problem that I need help with (thank you).

R = 0.1;
x0 = 0;
y0 = 0;
z0 = 0;
SubLengthX=0.001;
SubLengthY=0.001;
SubLengthZ=0.001;
z1=Re[z0-Sqrt[R^2-(x-x0)^2-(y-y0)^2]];
z2=Re[z0+Sqrt[R^2-(x-x0)^2-(y-y0)^2]];
y1=Re[y0-Sqrt[R^2-(x-x0)^2]];
y2=Re[y0+Sqrt[R^2-(x-x0)^2]];
x1=-R+x0;
x2=+R+x0;
Integrate3D[functionToIntegrate_,x1_,x2_,y1_,y2_,z1_,z2_,subLenX_,subLenY_,subLenZ_]:=Module[{volumeList,compiledSum,subVolume},
Print["Creating mesh..."];
volumeList=Partition[Flatten[ParallelTable[
Table[
Table[
{x+Mod[x2-x1,SubLengthX]/2,y+Mod[y2-y1,SubLengthY]/2,z+Mod[z2-z1,SubLengthZ]/2}
,Evaluate[{z,z1+SubLengthZ/2,z2-SubLengthZ/2,SubLengthZ}]]
,Evaluate[{y,y1+SubLengthY/2,y2-SubLengthY/2,SubLengthY}]]
,{x,x1+SubLengthX/2,x2-SubLengthX/2,SubLengthX}]],3];
Print["Compiling integrator..."];
compiledSum=Compile[{{subVol, _Real},{volumeList1,_Real,2}},
subVol ParallelSum[
(functionToIntegrate/.{x->volumeList1[[p,1]],y->volumeList1[[p,2]],z>volumeList1[[p,3]]})
,{p,Length[volumeList1]}]
,CompilationTarget->"C"];
Print["Integrating..."<>" Number of sub-volumes is "<>ToString[Length[volumeList]]];
subVolume=subLenX subLenY subLenZ;
compiledSum[subVolume,volumeList]
]

Thank you for your efforts.

closed as off-topic by Michael E2, xzczd, Bob Hanlon, bbgodfrey, Dr. belisariusApr 9 '15 at 13:29

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Michael E2, xzczd, Bob Hanlon, bbgodfrey, Dr. belisarius
If this question can be reworded to fit the rules in the help center, please edit the question.

• Why not to use built-in NIntegrate function? There are a lot of examples of its use in documentation. The domain of integratinon can be expressed by Ball function. – Ivan Apr 9 '15 at 10:41

For your example, I would simply do this:

R = 2.3;
{x0, y0, z0} = {1.2, 2.3, 3.4};
NIntegrate[1, {x, -R + x0, R + x0}, {y, y0 - Sqrt[R^2 - (x - x0)^2],
y0 + Sqrt[R^2 - (x - x0)^2]}, {z,
z0 - Sqrt[R^2 - (x - x0)^2 - (y - y0)^2],
z0 + Sqrt[R^2 - (x - x0)^2 - (y - y0)^2]}];

For a general function func = Function[{x,y,z},body] and a set of boundaries bound = {{x1,x2},{y1,y2},{z1,z2}}, I would do

integral3d[func_, bound_] :=

NIntegrate[
func[x, y, z], {x, bound[[1, 1]], bound[[1, 2]]}, {y, bound[[2, 1]],
bound[[2, 2]]}, {z, bound[[3, 1]], bound[[3, 2]]}]

The integral over a spherical region is easily performed by Mathematica even analytically. Assuming f=1 and for brevity putting the center of the sphere at the origin:

Timing@Integrate[1, {x, -r, +r}, {y, -Sqrt[r^2 - x^2], +Sqrt[
r^2 - x^2]}, {z, -Sqrt[r^2 - x^2 - y^2], +Sqrt[r^2 - x^2 - y^2]}]

(*{0.218401, 4 Pi r^3 / 3}*)

Assigning a numerical value to r:

r = 1;
Timing@Integrate[1, {x, -r, +r}, {y, -Sqrt[r^2 - x^2], +Sqrt[
r^2 - x^2]}, {z, -Sqrt[r^2 - x^2 - y^2], +Sqrt[r^2 - x^2 - y^2]}]

(*{0.171601, 4 Pi / 3}*)

Integrating numerically

Timing@NIntegrate[1, {x, -r, +r}, {y, -Sqrt[r^2 - x^2], +Sqrt[
r^2 - x^2]}, {z, -Sqrt[r^2 - x^2 - y^2], +Sqrt[r^2 - x^2 - y^2]}]

(*{0.015600, 4.18879}*)

Checking the result:

N[4 Pi / 3]

(*4.18879*)

PS Maybe I haven't understood the exact trouble you were facing, but if I have, the short answer is that Mathematica does have built-in functional for n-dimensional integration, i.e.

Integrate[f[x1,x2,...,xn],{x1,x10,x11},{x2,x20,x21},...,{xn,xn0,xn1}]

EDIT

For the spherical region not at origin specifying r,x0,y0,z0 with a list of replacement rules at the end:

NIntegrate[1, {x, -r + x0, +r + x0},
{y, y0 - Sqrt[r^2 - (x - x0)^2], y0 + Sqrt[r^2 - (x - x0)^2]},
{z, z0 - Sqrt[r^2 - (x - x0)^2 - (y - y0)^2],
z0 + Sqrt[r^2 - (x - x0)^2 - (y - y0)^2]}] /.
{x0 -> 1.2, y0 -> 0.3, z0 -> -0.6, r -> 3.}

Analytical integration also works, given the correct assumptions although the addition of extra undefined symbols makes Mathematica take much longer to compute:

Integrate[1, {x, -r + x0, +r + x0},
{y, y0 - Sqrt[r^2 - (x - x0)^2], y0 + Sqrt[r^2 - (x - x0)^2]},
{z, z0 - Sqrt[r^2 - (x - x0)^2 - (y - y0)^2],
z0 + Sqrt[r^2 - (x - x0)^2 - (y - y0)^2]},
Assumptions -> x0 \[Element] Reals && y0 \[Element] Reals && z0 \[Element] Reals && r > 0]

As a matter of taste shift the origin within the function being integrated rather than in the integration limits:

(N)Integrate[f[x-x0,y-y0,z-z0], {x, -r, +r},
{y, -Sqrt[r^2 - x^2], +Sqrt[r^2 - x^2]},
{z, -Sqrt[r^2 - x^2 - y^2], +Sqrt[r^2 - x^2 - y^2]}]

EDIT2 - regarding usage of spherical coordinates

Here's a simple example where despite its simplicity Mathematica just can't get around the integral in Cartesian coordinates. The following two definitions are equivalent:

f[x_, y_, z_] := Exp[-x^2 - y^2 - z^2]
f[r_] := Exp[-r^2]

Executing

Integrate[f[x, y, z], {x, -r, +r}, {y, -Sqrt[r^2 - x^2], +Sqrt[r^2 - x^2]},
{z, -Sqrt[r^2 - x^2 - y^2], +Sqrt[r^2 - x^2 - y^2]},
Assumptions -> r > 0]

returns

Sqrt[Pi]*Integrate[E^(-x^2 - y^2)*Erf[Sqrt[r^2 - x^2 - y^2]], {x, -r, r},
{y, -Sqrt[r^2 - x^2], Sqrt[r^2 - x^2]}, Assumptions -> r > 0]

i.e. Mathematica only takes the outermost integral. However, the equivalent integral in spherical coordinates

Integrate[r^2 Sin[th] f[r], {r, 0, r1}, {th, 0, Pi}, {phi, 0, 2 Pi},
Assumptions -> r > 0]

is evaluated fully and returns

(-2*Pi*r1)/E^r1^2 + Pi^(3/2)*Erf[r1]

With an arbitrary function defined in terms of Cartesian coordinates f[x,y,z] we can always plug f[r Sin[th] Cos[phi], r Sin[th] Sin[phi], r Cos[th] ] (with the r^2 Sin[th] multiplier, of course) into the integral. Or (with integration over a sphere not at the origin) like so:

(N)Integrate[f[r Sin[th] Cos[phi] - x0,
r Sin[th] Sin[phi] - y0,
r Cos[th] - z0] * r^2 Sin[th],
{r, 0, r1}, {th, 0, Pi}, {phi, 0, 2 Pi},
Assumptions -> r > 0]

Of course, the expanded form of

f[r Sin[th] Cos[phi] - x0,
r Sin[th] Sin[phi] - y0,
r Cos[th] - z0] * r^2 Sin[th]

will probably be horifically complicated, but chances are, that the engine will thanks us for the simple, constant integration boundaries.

• Sorry but this is a very bad answer. Notice that I'm not integrating at the origin necessarily, and you're welcome to integrate analytically over a sphere that's not at the origin. I'm -1ing this. – The Quantum Physicist Apr 9 '15 at 11:22
• As I said, I've omitted the not-at-origin for brevity. The substitution is trivial, for example: NIntegrate[ 1, {x, -r + x0, +r + x0}, {y, y0 - Sqrt[r^2 - (x - x0)^2], y0 + Sqrt[r^2 - (x - x0)^2]}, {z, z0 - Sqrt[r^2 - (x - x0)^2 - (y - y0)^2], z0 + Sqrt[r^2 - (x - x0)^2 - (y - y0)^2]}] /. {x0 -> 1.2, y0 -> 0.3, z0 -> -0.6, r -> 3.} Of course NIntegrate requires numerical values, hence the replacement rules at the end. – LLlAMnYP Apr 9 '15 at 11:36
• It's not like that, because I don't want to calculate just the volume. The volume thing is just for testing accuracy (as explained in the question). I want to integrate over a nonlinear function inside the sphere. – The Quantum Physicist Apr 9 '15 at 13:05
• The general expression for integrating over a spherical region is included in my (edited) answer. How well and quickly Mathematica will handle it depends on the function itself. The integral for getting the volume illustrates that Mathematica can integrate both analytically and with good numeric accuracy. On a case by case basis I would strongly recommend transforming the integral into spherical coordinates (this would make the job easier as the bounds would no longer be interdependent). – LLlAMnYP Apr 9 '15 at 13:18
• Thanks for editing your answer. Now my my question is answered by you and by the other guy. I removed the -1. About the spherical coordinates, again, it's a good solution ONLY if your sphere is centered in the origin. Otherwise, spherical coordinates become waaaaay more complicated. Again, try writing the general sphere equation in sperical coordinates and you'll see how impossible it's. – The Quantum Physicist Apr 9 '15 at 13:49