20

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?

Siong Thye Goh
  • 155,038
entrelac
  • 4,791
  • 7
    This summation is said to be separable, as it can be split in a pure horizontal pass and a pure vertical one, as the array of coefficients is the tensor product of the vectors. The split scheme can be seen as a numerical equivalent of Fubini's theorem. But, sorry, all of this is obvious. –  Jun 10 '15 at 14:21
  • So, you are taking a Cartesian product of the composite rule with itself, then. – J. M. ain't a mathematician May 26 '17 at 14:34
  • This may be an interesting observation if you use software like Maple which makes use of the built in matrix manipulation operations? – jim Aug 23 '17 at 10:22
  • Well, this goes a couple of years late, but could you explore this more? I'm currently trying to implement this but have made no advancements – Bidon Mar 23 '19 at 18:29

1 Answers1

1

The two–dimensional Simpson rule is nothing other than applying the familiar 1D composite Simpson operator in each coordinate in succession. Start by fixing $y$ and regarding

$$ G(y)\;=\;\int_a^b f(x,y)\,dx $$

as a function of $y$. Approximating this inner integral by the composite Simpson rule on $[a,b]$ gives

$$ G(y)\approx \frac{h_x}{3}\sum_{i=0}^m s_i\,f(x_i,y), $$

where the weights $s_i$ alternate like $1,4,2,4,\dots,2,4,1$. Treat $G(y)$ itself as the integrand in the outer integral

$$ \int_c^d G(y)\,dy, $$

apply the composite Simpson rule again, and get

$$ \int_c^d G(y)\,dy \;\approx\; \frac{h_y}{3}\sum_{j=0}^n t_j\,G(y_j) \;=\; \frac{h_x\,h_y}{9} \sum_{i=0}^m\sum_{j=0}^n(s_i\,t_j)\,f(x_i,y_j), $$

where $t_j$ are again the Simpson weights $1,4,2,\dots,4,1$ in the $y$-direction. Since each coefficient in the double sum factors as $W_{i,j}=s_i\,t_j$, the full weight matrix is simply the outer product of the one dimensional weight vectors.

Then the Fubini’s theorem to assert that

$$ \iint_{[a,b]\times[c,d]}f(x,y)\,dx\,dy =\int_c^d\Bigl[\int_a^b f(x,y)\,dx\Bigr]dy, $$

and then replace each 1D integral by its composite Simpson approximation. Because composite Simpson is exact for any polynomial in $x$ of degree up to three (at each fixed $y$), and likewise exact in $y$ for fixed $x$, the resulting 2D quadrature rule is exact for all bivariate polynomials of degree at most three in each variable separately. In other words, by thinking of Simpson’s rule as a linear functional in $x$ followed by one in $y$, one arrives naturally at the tensor product structure that you mentioned, and the “vertical” and “horizontal” coefficient arrays multiply to produce the familiar 2D weight matrix.

Suraj Gavhale
  • 200