How to integrate over a discrete 2D surface in MATLAB?

Question

I have a function z = f(x, y), where z is the value at point (x, y). How may I integrate z over the x-y plane in MATLAB?

By function above, I actually mean I have something similar to a hash table. That is, given a (x, y) pair, I can look up the table to find the corresponding z value.

The problem would be rather simple, if the points were uniformly distributed over x-y plane, in which case I can simply sum up all the z values, multiply it with the bottom area, and finally divide it by the number of points I have. However, the distribution is not uniform as shown below. So I am actually asking for the computation method that minimises the error.

Answer 1

The currently accepted answer will only work for gridded data. If your data is scattered you can use the following approach instead:

scatteredInterpolant + integral2:

f = scatteredInterpolant(x(:), y(:), z(:), 'linear');
int = integral2(@(x,y) f(x,y), xmin, xmax, ymin, ymax);

This defines the linear interpolant f of the data z(i) = f(x(i),y(i)) and uses it as an argument to integral2. Note that ymin and ymax, instead of doubles, can be function handles depending on x. So usually you will be integrating rectangles, but this could be used for integration regions a bit more complicated.

If your integration area is rather complicated or has holes, you should consider triangulating your data.

DIY using triangulation:

Let's say your integration area is given by the triangulation trep, which for example could be obtained by trep = delaunayTriangulation(x(:), y(:)). If you have your values z corresponding to z(i) = f(trep.Points(i,1), trep.Points(i,2)), you can use the following integration routine. It computes the exact integral of the linear interpolant. This is done by evaluating the areas of all the triangles and then using these areas as weights for the midpoint(mean)-value on each triangle.

function int = integrateTriangulation(trep, z)
P = trep.Points; T = trep.ConnectivityList;
d21 = P(T(:,2),:)-P(T(:,1),:);
d31 = P(T(:,3),:)-P(T(:,1),:);
areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
int = areas'*mean(z(T),2);

Answer 2

If you have a discrete dataset for which you have all the x and y values over which z is defined, then just obtain the Zdata matrix corresponding to those (x,y) pairs. Save this matrix, and then you can make it a continuous function using interp2:

function z_interp = fun(x,y)

    z_interp = interp2(Xdata,Ydata,Zdata,x,y);

end

Then you can use integral2 to find the integral:

q = integral2(@fun,xmin,xmax,ymin,ymax)

where @fun is your function handle that takes in two inputs.