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?

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.

• Let me see if I got it, basically the order remains the same for any dimension since every term we sum in any trapezoidal rule integral has to be evaluated N times and so on. One final question , since this is true we can conclude that this method has limitations in terms of accuracy since we can't get any better than $O(N^{-2})$ and the only way to get better precision is increasing $N$ right? May 20 '16 at 19:49
• That is basically correct, but let me add to the answer to make sure it's clear.
– RRL
May 20 '16 at 20:05
• Yes - you have to increase $N$ in each dimension to get a significant improvement in precision. When the number of dimensions gets to around $10$ this becomes intractable. Monte Carlo methods are far more efficient for many dimensions.
– RRL
May 20 '16 at 20:47
• @RRL Why does the fact that $N^d$ evaluations are required means that the accuracy is $\mathcal{O}(n^{-2/d})$? I understand that if your mesh width is $1/N$ for one dimensional that you need $N^d$ evaluations but I don't understand how this translates to the accuracy. Feb 14 '18 at 7:13