Error for Trapezoidal Rule in multi-variable integrals For one dimension integrals $\int_{a}^{b}f(x)dx $, we know the global truncation error goes like$\ \approx\mathcal{O}(h^2)$ where $h=\frac{b-a}{N}$ and N is the number of intervals. Also knowing how related $h$ and $N$ are we can rewrite the global error $\mathcal{O}(\frac{1}{N^2})$ if $N >> (b-a)$ (most of the times this is correct since we want high precision so we have a large $N$).
Now if we add one more integral ( for simplicity i only assume rectangular regions ) : $\int_{a}^{b}\int_{c}^{d}f(x,y)dydx$ and applying the Trapezoidal rule twice I'm not sure if we can say the global error goes like $\mathcal{O}(\frac{1}{N^4})$ ( also i assume the same number of intervals for $x$ and $y$).
Question: Is the global truncation error for a $n$ dimensional integral using the Trapezoidal rule $\mathcal{O}(\frac{1}{N^{2n}})$, assuming only "rectangular" regions of integration and every single integral has the same $N$  number of intervals?
 A: No.
Using product quadrature rules for multi-dimensional integrals suffers from the so-called curse of dimensionality.  An $O(N^{-2})$ accurate rule using N intervals in one-dimension is generally $O(N^{-2})$ accurate when applied as a product rule in $d$ dimensions. However there will be $M = N^d$ evaluations required. So the accuracy is $O(M^{-2/d}).$ 
Again, if $N$ is the number of intervals in each dimension, then $O(N^{-2})$ accuracy is generally attained assuming the function is well-behaved.   Viewed in terms of the total number of quadrature points, however, the order of accuracy gets much worse
To see how the order of accuracy with respect to $N$ carries through consider a simpler mid-point rule. 
If $f:[0,1] \rightarrow \mathbb{R}$ and the second-derivative is bounded -- $|f''(x)|\leq K$ -- then the error bound for the mid-point rule applied at the mid-points $\bar{x_i}$ of $N$ uniform invervals is
$$\left|\int_{0}^{1}f(x)\,dx - \frac1{N}\sum_{i=1}^{N}f\left(\bar{x_i}\right)\right|\leq \frac{K}{24N^2}.$$
In the two-dimensional case,  $f:[0,1]^2 \rightarrow \mathbb{R}$, where $f$ has bounded second-order partial derivatives we have the same order of accuracy with respect to N. Using a Taylor series approximation, we find that the first- and mixed second-order derivative terms vanish after integrating and the remaining error is
$$\left|\int_{0}^{1}\int_{0}^{1}f(x,y)\,dx\,dy - \frac1{N^2}\sum_{i=1}^{N}\sum_{j=1}^{N}f\left(\bar{x_i},\bar{y_i}\right)\right|\\=\left|\sum_{i=1}^{N}\sum_{j=1}^{N}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left[f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right]\,dx\,dy \right| \\ = N^2O(N^{-4})= O(N^{-2})$$
To provide more detail, using the triangle inequality, the global error bound is
$$E_N \leqslant \sum_{i=1}^{N}\sum_{j=1}^{N}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left|f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right|\,dx\,dy.$$
Assume that the second-order partial derivatives are bounded in absolute value by $K$ and the points are uniformly spaced with $x_i - x_{i-1} = 1/N$ and $y_j - y_{j-1} = 1/N.$
The local error bound on an individual rectangle is
$$E_{ij} = \\ \int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left|f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right| \, dx \, dy \\ \leqslant \frac{K}{2}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}(x-\bar{x_i})^2 \, dx \, dy +  \frac{K}{2}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}(y-\bar{y_i})^2 \, dx \, dy  \\ = \frac{K}{2}(y_j -y_{j-1})\int_{x_{i-1}}^{x_i}(x-\bar{x_i})^2 \, dx + \frac{K}{2}(x_i -x_{i-1})\int_{y_{j-1}}^{y_j}(y-\bar{y_i})^2 \, dy  \\ = \frac{K}{6N}\left.(x-\bar{x_i})^3\right|_{x_{i-1}}^{x_i} + \frac{K}{6N}\left.(y-\bar{y_j})^3\right|_{y_{j-1}}^{y_j}.$$
With $\bar{x}_i = (x_i + x_{i-1})/2$ and $\bar{y}_i = (y_j + y_{j-1})/2$ we have 
$$x_i - \bar{x}_i = -(x_{i-1} - \bar{x}_i) = (x_i - x_{i-1})/2 = 1/2N, \\ y_j - \bar{y}_j = -(y_{j-1} - \bar{y}_j) = (y_j - y_{j-1})/2 = 1/2N,$$
and
$$E_{ij} \leqslant  \frac{K}{6N}2\left(\frac{1}{2N}\right)^3 + \frac{K}{6N}2\left(\frac{1}{2N}\right)^3 = \frac{K}{12N^4}.$$
Summing over $N^2$ rectangles, the global error bound is 
$$E_N = N^2 \frac{K}{12N^4} =  \frac{K}{12N^2}$$
Hence, the local error is $O(N^{-4})$, but the global error is $O(N^{-2})$ as in the one-dimensional case.
