# Integration over a complicated domain

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$.

• Do you have the ability to determine for a given $x$ how many values the curve $C$ has? – Wolfgang Bangerth Apr 6 '15 at 15:14
• @WolfgangBangerth Yes, I do. At least empirically I can calculate the number of crossings. – Nick P Apr 6 '15 at 17:23
• If the function f is known you can reformulate the integration on the domain to an integration on the boundary with Stokes Theorem. Depending on an adequate parametrization for x=x(t), y=y(t) the integral of the bounding curve can be solved by quadrature techniques. – Bort Apr 7 '15 at 8:55
• @Bort: I think this is a bit misleading. Knowing $f$ does not guarantee that you can find a parametrization (easily). If you can find such a parameterization, however, your answer is definitely worth considering – bigge Apr 7 '15 at 19:04

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.

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$.

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.

• I've done some work on this approach, and you can get extremely high-accuracy results if you have some knowledge about the form of the bounding curve. The idea is that you can tell a quadrature routine about critical points where the function jumps (the curve C, in this case), and then use these to eliminate the discontinuity error, along with the inefficiencies. I have a project that addresses this specific problem at github.com/woodscn/BACL-MMS. It's still in prototypical shape, but I'd be interested to try it out on this problem. – woodscn Apr 9 '15 at 21:22
• Do you have something written down on the approach? – bigge Apr 10 '15 at 4:26
• Behind paywalls, unfortunately: researchgate.net/publication/… – woodscn Jun 1 '16 at 20:18