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