# Numerical integration of a three dimensional array

I need to integrate a scalar valued function $f\left(\boldsymbol{x}\right)$ where $\boldsymbol{x}$ is a three dimensional position vector; in other words:

$\int_{\Omega} f\left(\boldsymbol{x}\right) d^{3}\boldsymbol{x}$

The integration domain $\Omega$ is a a cube of side $L$.

The problem is that i don't have the analytical function but rather its values on a cartesian cubic grid $f \left(\boldsymbol{x_{\text{grid}}}\right)$. Specifically, in my program $f \left(\boldsymbol{x_{\text{grid}}}\right)$ is a three dimensional array.

I was thinking to write a program in which i use Gaussian quadrature with $n$ points inside each element and interpolating the value of the function in the integration points and sum up the result from each element. This is however quite an intensive work, therefore my question is:

Is there a straightforward or easier way to perform such an integration?

• Generate an InterpolatingFunction using Interpolation and then use NIntegrate. Alternatively, you could just Sum over f. Feb 24 '15 at 17:01
• Are your samples evenly spaced inside $\Omega$? Feb 24 '15 at 17:27
• By evenly spaced you mean equidistant one to each other? If so, yes Feb 24 '15 at 17:31
• @bbgodfrey I don't think I can just sum up f, isn't there a Volume involved in each summation? Feb 24 '15 at 17:40
• If the points are uniformly spaced, the volumes are the same for each point. So, multiply the total by the volume for one point; i.e., the volume of the cube divided by the number of points. Feb 24 '15 at 17:42

To elaborate on my Comment (and assuming uniform spacing of the data), consider the toy problem in 1-D:

f = Table[Sin[2 Pi ( i - .5)/10], {i, 10}]


Generating an InterpolatingFunction and then using NIntegrate yields:

g = Interpolation[f]
NIntegrate[g[x], {x, 1, 10}]
(* 3.3306690738754696*^-16 *)


Simply forming the Total yields the same result to machine precision.

Total[f]
(* 0. *)


This is not surprising, because using splines, etc cannot introduce more accurate information, because there is no more information. (Mathematically, the splines, etc integrate to one.)

On the other hand, if the points are not evenly spaced, then Interpolation plus NIntegrate is a straightforward (but not the only) approach.

• +1 the Interpolation approach is also useful for handling the edge conditions that arise depending on how your grid aligns with the region boundary. Feb 24 '15 at 17:53
• actually my grid is perfectly aligned with the edge boundaries, because it is defined on the boundaries (a cube side divided in $N$ segments) Feb 24 '15 at 17:58