Integration over a complicated domain

Question

Consider points $(x,y)$ in a domain $\mathcal{D}$ that is bounded above by a curve $C$ that is multivalued (i.e. I have points defining it that are of the form $(x_n,y_n)$, $n=1,2,..,N$), and looks like, for instance the figure shown below.

The remaining parts of the boundary $\partial \mathcal{D}$ occur at $x=0, x=L$ and $y=-h$.

How do I numerically compute general surface integrals of the form

$$\int_{\mathcal{D}} f(x,y) \ dx\ dy$$

for a function $f(x,y)$. For example, $f(x,y) = x^2+y^2$.

Answer 1

You could triangulate the domain (with say a Delaunay method), and then integrate over the triangles using standard integration techniques. This would be a finite element method, but just for integrating a function rather than solving a differential equation.

Answer 2

To give you a better picture, what I mean with my comment, consider

Stokes Theorem

$$\int_D db^{n-1}=\int_{\partial D}b^{n-1}$$

be $b^{n-1}$ an arbitrary 1-form $$b^{1}= a_1 dx + a_2 dy$$ leads to $$\int_D f(x,y)dxdy = \int_D\left(\frac{\partial a_2}{\partial x}-\frac{\partial a_1}{\partial y}\right)dxdy$$ Without loss of generality, $a_1=0$ one needs to solve $f(x,y)=\frac{\partial a_2}{\partial x}$ for $a_2$.

Now you need the parametrization for the boundary functions. For the first three bounding functions $$(x=0,y\in[-h,\eta(0)]),(x\in[0,L],y=-h),(x=L,y\in[-h,\eta(L)])$$ this is straight forward. For the bounding function I assume you know or can evaluate the function in $(x(t),y(t)) \wedge t\in[0,1]$. So I suggest using two Gaussian quadratures. First from $[lb,peak]$ and second from $[peak, rb]$.

Given you can obtain information on the peak and integrate $f$, this method will give you very high accuracy and very good convergence.

Supplementing Edit

As bigge pointed out, knowing $f$ does not guarantee to find an parametrization easily. Of course this is true. But even if one is not able to solve $f(x,y)=\frac{\partial a_2}{\partial x}$ analytically, one can still use this approach. Assuming Gaussian quadrature techniques one can read it as

$$a_2(x,y)=\int_{x_0}^x f(t,y)dt =\int_{-1}^1 f(r,y) \frac{dt}{dr} dr\approx \sum_k w_k f(r_k,y) \frac{dt}{dr}$$ With e.g. $x_0=0$ chosen arbitrarily. Now the scaling factor $\frac{dt}{dr}$ depends on the integration range $x$. Since we still integrate line integrals along the boundary, $x$ is not arbitrary but specific set of quadrature points for evaluating $\int_{\partial D} a_2(x,y) dy$.
So the final term is a double sum for each boundary integral as expected for an integral over a 2-D domain. In that form it is simply a very sophisticated 2D Simpson's rule.

In my opinion, it is usually more challenging to find a good parametrization of the bounding curves than evaluation of $f$.

Answer 3

Depending on your accuracy requirements and the effort you are willing to put into the method, you can do different things, especially since you already have an implicit description of your surface. Here's the idea: Using the example of the integral over a circle, you can define a function $$\varphi(x,y) = 1 - (x^2 + y^2)$$ that is positive inside of your domain of integration (in this case, a circle with radius 1) and negative anywhere else. $\varphi(x,y)$ is then called a level set function and there's tons of papers dealing with the integration over the enclosed area/volume $(\varphi(x,y) > 0)$ and over the corresponding domain boundary $(\varphi(x,y) = 0)$. The corresponding integrals can simply be rewritten as $$\int_D H(\varphi(x,y)) \,dV$$ and $$\int_D \delta(\varphi(x,y)) \,dV,$$ respectively, where $H$ is the Heaviside function and $\delta$ is the Dirac delta distribution.

These integrals cannot be computed directly, but they can be evaluated using different methods which strongly differ in terms of complexity and accuracy. One approach has already been described by Bill, but, if your accuracy requirements are relatively low, you can get away with a few very simple options:

• Use simple adaptive quadrature techniques (e.g., such the ones implemented in Matlab)

• Replace $H$ and $\delta$ with smoothed approximations (see Wikipedia for a few examples) and evaluate them using classical Gauss quadrature

• Monte-Carlo

If your accuracy requirements are high, you will have to invest a little more into the topic, since the above-mentioned approaches tend to become exceedingly inefficient when you increase the resolution (they just converge very slowly). In this case, it's not easy giving a simple answer without knowing a little more about your problem.