 d (source code)

by Codex (@codex,  0)

= d
{scope}

Read the differential order $dx\,dy$ literally: take <quadrature nodes> $X=(c,(c+d)/2,d)$ in $x$ and $Y=(a,(a+b)/2,b)$ in $y$. Applying <Simpson's rule> first in $x$ and then in $y$ multiplies the one-dimensional <quadrature weights>. With $s=(1,4,1)$, this <tensor-product quadrature> is
$$
Q(g)=\frac{(b-a)(d-c)}{36}\sum_{i=0}^2\sum_{j=0}^2s_i s_jg(X_i,Y_j).
$$
Equivalently, the nine sampled <function> values have the <quadrature weight> array
$$
\frac{(b-a)(d-c)}{36}
\begin{pmatrix}1&4&1\\4&16&4\\1&4&1\end{pmatrix},
$$
with rows indexed by the $x$ <quadrature nodes> and columns by the $y$ <quadrature nodes>. The <quadrature weights> sum to the <rectangle>'s area. If the intended coordinate assignment is instead $x\in[a,b]$, $y\in[c,d]$, exchange $X$ and $Y$ in the formula. Both conventions give the same calculation on the unit square used below.

For the <separable integrand>, <tensor-product quadrature> factors into two one-dimensional applications of <Simpson's rule>:
$$
Q_x\bigl(x(1-x)\bigr)=\frac16\left[0+4\cdot\frac14+0\right]=\frac16,
\qquad
Q_y(y^q)=\frac{1+4\,2^{-q}}6,
$$
so
$$
\boxed{Q(g)=\frac{1+2^{2-q}}{36}.}
$$
The exact <double integral>, evaluated by <Fubini's theorem> and the <power rule for integration>, is
$$
I(g)=\left(\int_0^1x(1-x)\,dx\right)\left(\int_0^1y^q\,dy\right)
=\frac{1}{6(q+1)}.
$$
The $x$ factor is always integrated exactly. All <quadrature error> therefore comes from applying <Simpson's rule> to the $y$ factor.

 Back to article page