# Numerical Integration over 2D Regions with Discontinuous Functions

I've run into a tricky problem, and I haven't really thought up a good solution. I have to compute many integrals of the form

$$\iint\limits_{B((x_k,y_k),\epsilon))} f(x,y) dy dx$$ where $B((x_k,y_k),\epsilon)$ is a ball about the point $(x_k,y_k)$ of radius $\epsilon$. The function $f$ takes values in $[0,1]^2$, and is zero outside. For $(x_k,y_k)$ sufficiently in the interior of the unit square, this isn't so bad and I can apply a linear change of variables by choosing $x' = x-x_k$ and $y' = y-y_k$, then going to polar coordinates. Once I'm in polar coordinates, I can call a routine to do Gauss-quadrature by creating quadrature points on $[0,\epsilon]$ and $[0,2\pi]$, and then taking the Cartesian product of the 1D nodes and products of the weights to get a 2D quadrature scheme over the region in polar coordinates. This works fairly well. If anyone has suggestions for alternatives, I am more than interested in hearing of any other ideas.

Here is the trouble I'm having: If $f$ is continuous and goes to zero smoothly along the boundary of the square, then the above scheme works fairly well (recall that I have declare $f$ to be zero outside of the unit square). My problem is the case when $f$ does NOT go to zero along the boundary the square. Then, for $(x_k,y_k)$ near the boundary, we have an integration region of $B((x_k,y_k),\epsilon)$, but for certain points $(x_k,y_k)$ the integration region will contain area outside of the unit square. There is a discontinuity along the boundary from the unit square $[0,1]^2$ to the outside, and suddenly my polar coordinate Gaussian quadrature starts performing very poorly. Since $f$ is zero outside of $[0,1]^2$, we can simply say our integration region is just $B((x_k,y_k),\epsilon) \cap [0,1]^2$, but then the question is how to easily integrate over this region? I have many integrals to compute and this is a bit troubling. I can't just do my polar coordinate change + product of Gauss nodes anymore, since my region is the intersection of the square with the circle.

A further difficulty is added when $f$ is discontinuous inside of $[0,1]^2$. Then I'm in a lot of trouble; even for $(x_k,y_k)$ sufficiently in the interior such that we don't have to worry about the discontinuity along the edge of the unit square, the polar-coordinate plus Gaussian 2D quadrature performs poorly because of the discontinuity. (For an example, consider a function $f(x,y)$ which is 1 for $0 \leq x \leq \frac{1}{2}$ and for $0 \leq y \leq 1$, and zero everywhere else. For $\epsilon = \frac{1}{8}$, the balls $B((x_k,y_k),\epsilon)$ often go over the discontinuity and the quadrature I'm doing performs poorly).

I am open to any suggestions or references. One thing I considered was possibly "meshing" the integration region since it is a circle, or a circle intersected with the square. I could "polygonize" the circle intersected with the square and triangulate this region. and and do some form of quadrature on each triangle. In particular, I'm looking for something that is implementable, and ideally not particularly slow. I attempted to use MATLAB's integral2 command and define my functions and regions, and it performed very slowly, and often returned a "failed to pass the global error test" message.

Please let me know if anything is unclear or confusing, and I will try to clarify. Thanks for the help!

• If your discontinuity is a line then your circle is cut into a region that can be broken into a sector plus a triangle. Perhaps you can integrate separately on the sector by using your polar version and on the triangle by using some Cartesian quadrature. Jul 25, 2013 at 4:05
• Can you give an example of the function $f(x,y)$ under consideration? Jul 25, 2013 at 4:10
• Maesumi: I considered that possibility of splitting along a line (for a step function example I listed, that's exactly the case). I think you're exactly right, and I should probably implement the sector with Gaussian quadrature I was using + a quadrature method over the triangle. I was partly concerned because I have to do this for many $(x_k,y_k)$ pairs. I need to figure out a good implementation, and maybe I can use a GPU or something to speed it up user64494: See my third paragraph for an example of a step function and an epsilon. Jul 25, 2013 at 4:35
• Do the balls have to be balls? What about using square boxes for integration that can be aligned and matched with the boundary? Jul 25, 2013 at 9:40
• BTW, your polar weights should include the radius as a factor. Jul 25, 2013 at 10:08

My suggestion is to use Maple to this end. For example, the Maple code $$f := unapply(piecewise(x >= 1/2 \,and\, x <= 1\, and \,y >= 1/2 \,and\, y <= 1, 1, 0), x, y):$$ $$a:=0.45:b:=0.5:$$ $$VectorCalculus:-int(f(x, y), [x, y] = Circle(<a, b>, .125), numeric, epsilon = 10^{-3})$$ produces the answer $0.006196086871$ with relative error $10^{-3}$ in 8 min on my comp, where another hard problem is simultaneously solving.
• Thanks for the suggestion. Unfortunately, I do not have Maple, and I also need to do several thousand of these integrals (each integral over a different ball $B((x_k,y_k),\epsilon)$ where there are thousands of points $(x_k,y_k)$, so 8 minutes per integral is a bit lengthy. Jul 28, 2013 at 19:30