# Numerical integration of a multivariate function whose expression is unknown

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.

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.

• Nice ! I should plot some examples of $l(x,y)$ to see what it looks like. – bela83 May 5 '15 at 5:54
• @bela83 I added a plot to show how $l(x,y)$ is constructed. If necessary, you can also derive other higher-order quadrature rules by using higher-order interpolating polynomials. – Kirill May 5 '15 at 6:09
• Can you elaborate on the plot ? It seems like the values at the corners are always $1$ (or close to), so I don't understand why the interpolant is not constant. Also I wanted to edit your post to correct the small typo : "the integral of $l(x,y)$ will be" – bela83 May 5 '15 at 7:16
• @bela83 See edit. I don't think I understand what you mean by "why the interpolant is not constant". – Kirill May 5 '15 at 7:25
• What values did you choose for $f(0,0), f(0,1), f(1,0), f(1,1)$ ? Or maybe you plotted the elementary interpolants ? OK I get it now. – bela83 May 5 '15 at 7:30

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.

• Thanks. 1 and 2 are the same, right ? – bela83 May 5 '15 at 5:56
• Conceptually, they're different. But I hadn't thought to check if they gave different results. After checking a few different thoughts, they all turn out to give the same results (which is the equation in Aurelius' answer). – Brendan May 5 '15 at 20:36

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.

• No I can't sample $f$ at arbitrary locations. – bela83 May 4 '15 at 16:06
• Then your options are to do the simple summation I wrote above, which will give accuracy comparable to trapezoidal rule, or to a multidimensional lagrange polynomial if you need higher accuracy. – Aurelius May 4 '15 at 16:11