Numerical integration of a multivariate function whose expression is unknown

Question

This question is related to the previous one, but I believe the extension to several variables is a problem in itself.

I have a collection of points $(x_j,y_k,f(x_j,y_k))$ for $x_j=\frac{j}{n}$ and $y_k=\frac{k}{n}$ a grid over $[0,1]^2$. These values come from a PDE solver for $f$. I need to compute the integral of $f$ on $[0,1]^2$. How could I make it ?

In other words, what is the natural extension of the trapezoidal rule for higher dimensions ? As I see it, the problem is that, while there is a unique segment interpolating two points, there can be no plane interpolating four points.

Answer 1

There is a simple extension of trapezoid rule to 2d. In 1d, if you have the interval $[0,1]$, the rule takes the linear interpolant $$ l(x) = f(0)(1-x) + f(1)x, $$ and integrates it to get $$ \int_0^1 l(x) = \tfrac12 f(0) + \tfrac12 f(1). $$ To integrate over a large interval, you then sum up the contributions of the individual small intervals to get the usual trapezoid formula: $$ h \big(\tfrac12 f(x_0) + f(x_1) + \cdots + f(x_{n-1}) + \tfrac12 f(x_n) \big). $$

In 2d, on the square $[0,1]^2$, you can interpolate the function linearly on the bottom edge $y=0$, and separately on the top edge $y=1$, and then interpolate linearly between the two linear interpolants: $$ l_0(x) = f(0,0)(1-x)+f(1,0)x, \qquad l_1(x) = f(0,1)(1-x)+f(1,1)x, $$ $$ l(x,y) = l_0(x)(1-y) + l_1(x)y. $$ The resulting interpolant is a linear combination of elementary interpolants on the square $[0,1]^2$, with coefficients being the function values, and the elementary interpolants look like this:

where the elementary interpolants are the coefficients of the function values $f(x,y)$ in the interpolant constructed from them: $$ l(x,y) = f(0,0)(1-x)(1-y) + f(1,0)x(1-y) + f(0,1) (1-x)y + f(1,1)xy, $$ so of the four functions $xy$, $x(1-y)$, $(1-x)y$ and $(1-x)(1-y)$.

The integral of $l(x,y)$ will be $$ \int_{[0,1]^2}l(x,y)\,dx\,dy = \tfrac14\big( f(0,0)+f(0,1)+f(1,0)+f(1,1) \big). $$

Summing up these contributions over small intervals would then give $$ \frac{h^2}{4}\sum_{j,k} w_{j,k} f(x_j, y_k) $$ where $w_{j,k}$ is an integer that counts how many squares the point $(x_j,y_k)$ belongs to: one in the corners, two on the edges, four in the middle.

Answer 2

There's no reason to think that you'll find a plane that intersects all four points in general, but you can try various other approximations. A few quick ideas:

1. Some sort of midpoint rule where you assume that the average value of $f $ across the region $[j/n,(j+1)/n] \times [k/n,(k+1)/n]$ is simply the average of the four corner values.

2. Instead of assuming a plane through all four corners, assume a linear interpolation in $x$ along the slices $y = k/n$ and $(k+1)/n$, then a linear interpolation in $y$ between those two linear fits; this results in an equation of the form $a + bx + cy + dxy$.

Some Googling for things like "2D numerical integration" will turn up a variety of methods. It's all variations on the idea of estimating the value of $f$ between your data points and then integrating between the data points. Which method you choose will depend on how accurate you want your results to be, and how much work you want to do to compute the integral.

Answer 3

If you have a fixed grid, the simplest integration you can do would be to sum up the average values multiplied by the area they're on, like:

$$ \int_0^1 \int_0^1 f(x,y) dx dy \approx \sum_{i,j}^{NI-1,NJ-1} 1/4 (f_{i,j}+f_{i+1,j}+f_{i,j+1}+f_{i+1,j+1})\Delta x \Delta y$$

If you're just doing some simple post-processing, this may be sufficiently accurate.

In general if you can sample $f(x,y)$ at arbitrary locations, you'd want to look at doing multi-dimensional quadrature. There is free code for this in the GNU scientific library and others. Higher-dimensional quadrature is done by tensor products of the 1D case.

As I see it, the problem is that, while there is a unique segment interpolating two points, there can be no plane interpolating four points.

True, although you can always construct a unique Lagrange polynomial which you can then integrate analytically. Any of this is available in matlab.