Double integration over a polygon in Matlab

Question

I am given a function @f(x,y) and I want to evaluate the integral of this function over a certain convex polygon in MATLAB. The polygon is not necessarily a rectangle and that's why I can't use MATLAB's function "dblquad". The polygon I have is given by a set of vertices represented by the vectors X and Y, i.e. the vertices are (X(1),Y(1)),....,(X(n),Y(n)). Is there any function or method that I can use?

Answer 1

The trick is to use tools to integrate inside the region of interest. I've written a few tools for integration in a triangulated domain.

% Define a function to integrate.
% This function takes an nx2 array, where each row
% contains a single point to evaluate the kernel at.
% This computes x^2 + y^2 at each point.
fun = @(xy) sum(xy.^2,2);

% define the domain as a triangulated polygon
% this tool uses ear clipping to do so.
sc = poly2tri([1 4 3 1],[1 3 5 4]);

% Gauss-Legendre integration over the 2-d domain
[integ,fev]= quadgsc(fun,sc,2)
integ =
       113.166666666667
fev =
     8

% the triangulated polygon...
plotsc(sc,'facecolor','none','markerfacecolor','r')
axis equal
grid on

We can visualize the function itself, as a mapping z(x,y) over that polygonal domain. When a range field is supplied, the simplicial complex turns into a 2-1 mapping from the 2-d (x,y) domain.

sc2 = refinesc(sc,'max',.5);
sc2.range = fun(sc2.domain);
plotsc(sc2,'markerfacecolor','r')
grid on
view(17,12)

This is a simple polynomial function over the domain of interest, so the default low order Gaussian integration was adequate. The scheme used is a Gauss-Legendre one in a tensor product form over a triangle, not truly optimal, but viable. The problem with Gaussian quadrature, is it is not adaptive. It computes an estimate, based on implicit approximation by polynomials over a finite set of points.

The above estimate used 8 function evals to compute that estimate. Since the kernel is a low order polynomial, it should do perfectly. The problem is, you need to know if it is a correct solution. This is the problem with a Gaussian quadrature, there is no simple way to know if the answer is correct, except for resolving the problem with a higher order scheme until it seems to converge.

See that with 1 point per triangle at the barycenter, we get the wrong answer, but the higher order estimates all agree.

[integ,fev]= quadgsc(fun,sc,1)
integ =
       107.777777777778
fev =
     2

[integ,fev]= quadgsc(fun,sc,3)
integ =
       113.166666666667

fev =
    18

[integ,fev]= quadgsc(fun,sc,4)
integ =
       113.166666666667
fev =
    32

After writing quadgsc, I had to try an adaptive solver, that works in the same way as the other quad tools do in MATLAB. This does an adaptive refinement of the triangulation, looking for triangles where the solution is not stable. The problem is, I never did finish writing these tools to my satisfaction. There are many different methods one can employ for the cubature problem over a triangulated domain. quadrsc does a low order solution, then refines it, uses a Richardson extrapolation, then compares the results. For any triangles where the difference is too large, it refines them further until it converges.

For example,

[integ,fev]= quadrsc(fun,sc)
integ =
          113.166666666667

fev =
    16

So this works. The problem shows up on more complex kernels, where the issue becomes to know when to stop the refinement, and to do so before one has used up too many function evaluations. I never did get that fully working to my satisfaction, so I never posted these tools. I can send the toolbox to those who send me direct mail. The zip file is about 2.4 MB. One day I'll get around to finishing those tools, I hope...