Simpson's Rule for double integrals: $$\int_a^b\int_c^df(x,y) \,dx \,dy$$ is given by $$S_{mn}=\frac{(b-a)(d-c)}{9mn} \sum_{i,j=0,0}^{m,n} W_{i+1,j+1} f(x_i,y_j) $$ where: $$W= \begin{pmatrix} 1&4&2&4& \ldots&4&1\\ 4&16&8&16&\ldots&16&4\\ 2&8&4&8&\ldots&8&2\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\\ 1&4&2&4&\ldots&4&1\\ \end{pmatrix}$$ and where $m$ and $n$ are the subdivisions of the domains along $x$ and $y$ respectively. They are obviously both even.
The proof usually involves applying Simpson's Rule to the first integral $\int_a^bf(x,y)dx=g(x,y)$, then applying it to the second $\int_c^dg(x,y)dy$ thus obtaining the coefficients (or weights) $W$ of $f(x,y)$.
I've thought of an alternative: first we find the area under $f(x_i,y_j)$ where $i$ is a fixed value but $y_j$ goes from $c$ to $d$ using Simpson's Composite Rule and fill the coefficients in a vertical matrix $$ \begin{pmatrix} 1\\ 4\\ 2\\ 4\\ 2\\ \vdots \\ 4\\ 1\\ \end{pmatrix} $$
Similarly for a fixed $j$: $$\begin{pmatrix} 1&4&2&4&2& \ldots &4&1\end{pmatrix}$$
By multiplying the vertical matrix by the horizontal one we get $W$. I thought of the integral (surface) as a linear combination of the $f(x_i,y_j)$'s, justifying the product of the two matrices as a composition of both iterations of Simpson's rule for each variable. Does any of this make sense? Is there a more formal proof of these statements?
Also, can you recommend a textbook or some online resource where I can find more approximations like this one?
