error estimates for multi-dimensional Riemann sums Suppose that $f$ is a continuous function of bounded variation from $R^2$ to $R$ that's negative outside of some bounded set, 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$?  ($O(1/n)$ or $O(1/n^2)$ or what?)
Also how can this basic bound be improved if one knows more about $f$, e.g. that it is smooth or concave?
I'm restricting the question to functions on $R^2$ for definiteness, but I'd like to know the more general situation for $R^n$.
 A: With the hypotheses given, one can't do better than $O(1/n)$ decay.  Consider for instance the function $\frac{1}{n} \cos^2(2\pi n x_1)$ smoothly localised to a ball for some large $n$.  This has a total variation norm of $O(1)$, but for this specific value of $n$, the Riemann sum will be off by $O(1/n)$.
Of course, this function depends on $n$.  For an $n$-independent example, one could then consider the Weierstrass type function $\sum_{n=1}^\infty \frac{1}{j^2 n_j} \cos^2(2\pi n_j x_1)$ smoothly localised to the unit ball, where $n_j$ goes rapidly to infinity.  This is still continuous and of bounded variation, but now the Riemann sum will be off by about $O(1/j^2 n_j)$ at scale $1/n_j$.
In dimensions $d$ greater than 1, the situation is much worse; one can't do much better than $O(1)$, basically because of the failure of the Sobolev embedding $W^{1,1} \subset L^\infty$ in higher dimensions.  For instance, one can consider a function $f$ that consists of a bump function of height 1 localised to a ball of radius $O( n^{-d/(d-1)} )$ at each lattice point on $\frac{1}{n} {\bf Z}^d \cap B(0,1)$.  This has total variation norm $O(1)$ and is bounded by $O(1)$, but the Riemann sum is off by $O(1)$.  By superimposing several such examples together as in the Weierstrass type example we can then construct an $n$-independent function of bounded variation and continuous of compact support whose Riemann sum error decays as slowly as one pleases.
Once one does have enough regularity (in, say, a Sobolev class) to control local $L^\infty$ oscillation, then one can estimate the error term in the Riemann sum by partitioning space into cubes, using some sort of local Sobolev inequality on each cube, and summing up.  This for instance gives an $O(1/n)$ error term in the one-dimensional bounded variation case.
One can also analyse Riemann sums by Littlewood-Paley theory.  Functions whose Fourier transform is supported on frequencies much smaller than $n$ have excellent agreement between the integrals and their Riemann sums (particularly if one uses quadrature to improve the accuracy of the latter), and functions whose Fourier transform are supported on frequencies much larger than $n$ have a negligible integral.  So the error term is basically the same thing as the Riemann sum of the high-frequency component of the function $f$.
Concavity should be very helpful, ruling out the oscillatory counterexamples mentioned above and giving some new bounds on first and second derivatives of $f$ that can be plugged into the local Sobolev inequality method, but I don't immediately see what the best bounds would be with this hypothesis.
A: If $f$ is of bounded variation, then there are bounds given by a (multi-dimensional generalization of a) theorem of Koksma. A reference is Kuipers and Niederreiter, Uniform Distribution Of Sequences. 
EDIT: Here are a couple of results from that book. 
Theorem 5.5: Koksma-Hlawka Inequality. Let $f(x)$ be of bounded variation on $[0,1]^k$ in the sense of Hardy and Krause. Let $\omega$ be the finite sequence of points 
$${\bf x}_1,\dots,{\bf x}_N$$    
in $[0,1]^k$, and let $\omega_{j_1,m\dots,j_p}$ denote the projection of the sequence $\omega$ on the $k-p$-dimensional face of $[0,1]^k$ defined by $x^{(j_1)}=\cdots=x^{(j_p)}=1$. Then we have 
$$
\left|{1\over N}\sum_{n=1}^Nf({\bf x_{\it n}})-\int_{[0,1]^k}f({\bf x})d{\bf x}\right|\le\sum_{p=1}^k\sum_{1,\dots,k;p}^*D_N^*(\omega_{p+1,\dots,k})V^{(p)}(f(\dots,1,\dots,1))
$$ 
where $V^{(p)}(f(\dots,1,\dots,1))$ denotes the $p$-dimensional variation of $f(x^{(1)},\dots,x^{(p)},1,\dots,1)$ on $[0,1]^p$ in the sense of Vitali and where the term of the sum corresponding to $p=k$ is understood to be $D_N^*(\omega)V^{(k)}(f)$. 
Here $D$ is a discrepancy, probably very simple to calculate for the situation at hand, but I'm not up to typing it out. I'm going to bail on typing out Theorem 5.6, too; it applies when $f$ has certain continuous partial derivatives, and replaces the variation $V$ with an integral of the absolute value of said derivatives. 
