# 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. Commented Feb 24, 2015 at 17:01
• Are your samples evenly spaced inside $\Omega$? Commented Feb 24, 2015 at 17:27
• By evenly spaced you mean equidistant one to each other? If so, yes Commented Feb 24, 2015 at 17:31
• @bbgodfrey I don't think I can just sum up f, isn't there a Volume involved in each summation? Commented Feb 24, 2015 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. Commented Feb 24, 2015 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. Commented Feb 24, 2015 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) Commented Feb 24, 2015 at 17:58