Suppose that $f$ is a twice-differentiable concave function from $R^2$ to $R$ that's negative outside of some bounded set (e.g. $f(x,y)=1-x^2-y^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 http://mathoverflow.net/questions/72891/error-estimates-for-multi-dimensional-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 non-negative.
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 $N-x$ on $y>0, -y+u\leq x\leq y-u$ 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).