Derivation of 2D Trapezoid Rule According to this reference: http://mathfaculty.fullerton.edu/mathews/n2003/SimpsonsRule2DMod.html the trapezoid rule in two dimensions in a domain  $x\in[a,b]$, $y\in[c,d]$ is:
\begin{align}
\iint\limits_A f(x, y)\ dA &= \int\limits_a^b\int\limits_c^d f(x, y)\ dy\ dx\\
                           &={ \Delta x \Delta y \over 4}\bigg[f(a, c) + f(b, c) + f(a, d) + f(b, d) \\
&+ 2\sum_i f(x_i, c) + 2\sum_i f(x_i, d) + 2\sum_j f(a, y_j) + 2\sum_j f(b, y_j)\\ &+ 4\sum_j \big(\sum_i f(x_i, y_j)\big)\bigg]
\end{align}
Considering that the trapezoid rule in 1D is $\int f(x) dx = \sum_i {f(x_i) + f(x_{i+1})\over 2}\Delta x_i$ I tried to derive the 2D trapezoid rule as follows:
\begin{align}
\iint\limits_A f(x,y)\ dA &= {\Delta x \over 2}\int dy \sum_i f(x_i,y) + f(x_{i+1}, y)\\ &={\Delta x \over 2} \sum_i \int dy\ (f(x_i, y) + f(x_{i+1}, y))\\ &={\Delta x \over 2} \sum_i \bigg[\int dy \ f(x_i, y) + \int dy\ f(x_{i+1}, y)\bigg]\\ &= {\Delta x \Delta y \over 4}\sum_i\bigg[\sum_j f(x_i, y_j) + f(x_i, y_{j+1}) + \sum_j f(x_{i+1}, y_j) + f(x_{i+1}, y_{j+1})\bigg]
\end{align}
but this is significantly different from the result shown above so I'm looking for help to derive the 2D trapezoid rule. Thanks.
 A: The 2D rule you found on the other Web page is for a uniform grid
(equal spacing between all the $x$ coordinates in one list,
equal spacing between all the $y$ coordinates in the other list).
With a uniform grid, the 1D rule becomes
\begin{align}
\int f(x) dx &= \sum_{i=0}^{n-1} \frac{f(x_i) + f(x_{i+1})}{2}\Delta x_i \\
&= \frac{\Delta x}{2}\sum_{i=0}^{n-1} \left( f(x_i) + f(x_{i+1}) \right)
&& \mbox{all $\Delta x_i$ are equal} \\
&= \frac{\Delta x}{2} \left(\sum_{i=0}^{n-1} f(x_i)
       + \sum_{i=0}^{n-1} f(x_{i+1})\right)  && \mbox{addition is commutative} \\
&= \frac{\Delta x}{2} \left(\sum_{i=0}^{n-1} f(x_i)
       + \sum_{i=1}^{n} f(x_i)\right) && \mbox{renumber the index $i$} \\
&= \frac{\Delta x}{2} \left(f(x_0) + f(x_n) + 2\sum_{i=1}^{n-1} f(x_i)\right)
&& \mbox{combine equal terms.}
\end{align}
The last formula is the usual form of the uniform-grid formula,
except that people usually list the $f(x_n)$ term last inside the parentheses.
I arranged things inside the parentheses so that it is clearer that $f(x_n)$ is not part of the "$\sum$" sum.
Note that this is one of those cases where shorthands like "$\sum_i$" are likely to cause errors, because with the trapezoid rule it really does matter exactly where you start adding up the terms and where you stop.
You need to pay careful attention to the first and last index of each sum in order to make sense of the derivation above.
If you are not convinced by my formal manipulation of sums, try working out an actual trapezoid approximation with some number of terms, for example $x_1$ through $x_5,$ and see how it works.
The cited 2D rule is what you should get if you follow steps similar to the ones you attempted, but using the usual uniform-grid 1D rule in the form I derived above.
