Suppose that $f$ is a twicedifferentiable concave function from $R^2$ to $R$ that's negative outside of some bounded set (e.g. $f(x,y)=1x^2y^2$) and let $F=$max$(f,0)$. Let $S_n$ be the Riemann sum for the integral of $F$ over $R^2$ obtained by summing the values of $F$ at all points in the lattice $(Z/n)^2$ and dividing by $n^2$. What sort of bounds can be given for the difference between $S_n$ and the integral of $F$ over $R^2$? Is it $O(1/n)$ or $O(1/n^2)$ or what? This is a more focussed version of the question error estimates for multidimensional Riemann sums .

It looks like the error is in $O(1/n^2)$, with a precise and optimal bound $C/n^2$ if you have a fixed bound on (1) the second derivative of the function (2) the radius of the region where it is nonnegative. As the question is stated there are two sources for the error term:
(Note that a complete argument has to be more precise because the function $f$ could have zero derivative at the points where it vanishes, then the number of boundary squares is $O(n^2)$ but I think the result does not change). To check that this estimate is optimal you can think of a function which is invariant under a rotation of angle $\pi/2$ and equal to say $Nx$ on $y>0, y+u\leq x\leq yu$ for some small $u>0$. Then the first error term can be made smaller than the second, while the second "boundary" error term is indeed of the order of $1/n^2$ (the boundary errors all sum up). 

