Question 3 of 7: Gauss quadrature over a bilinearly mapped rectangle

Question 3: Gauss quadrature over a bilinearly mapped rectangle (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. The integrand $f(x,y) = x(x^{2}+y)$, the rectangle $2 \le x \le 6$, $1 \le y \le 7$, the four bilinear shape functions above, and the exact value $g = 2304$ for comparison.

Find. The Jacobian of the isoparametric map, a Gauss–Legendre evaluation of $g$ presented in vector–matrix form, and an explanation of how the numerical result compares with the exact one.

Approach. Map the rectangle onto the parent square with the given bilinear functions, evaluate the Jacobian (which is constant for a rectangle), choose the lowest rule that is exact for the mapped polynomial degree, and sum the weighted station values as a matrix triple product.

1. Map the geometry. With the corner coordinates ordered to match the given shape functions — node 1 at $(2,1)$, node 2 at $(6,1)$, node 3 at $(2,7)$, node 4 at $(6,7)$ — the isoparametric map is $$x(\xi,\eta) = \sum_{i=1}^{4} N_i x_i = 4 + 2\xi,\qquad y(\xi,\eta) = \sum_{i=1}^{4} N_i y_i = 4 + 3\eta$$ Every bilinear cross-term cancels because opposite edges are parallel: the map of a rectangle is affine.
2. Form the Jacobian. Differentiating the map, $$[J] = \begin{bmatrix} \dfrac{\partial x}{\partial \xi} & \dfrac{\partial y}{\partial \xi}\\[6pt] \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \eta}\end{bmatrix} = \begin{bmatrix} 2 & 0\\ 0 & 3\end{bmatrix},\qquad |J| = \det[J] = 6$$ Because $[J]$ contains no $\xi$ or $\eta$, the area element $dx\,dy = |J|\,d\xi\,d\eta = 6\,d\xi\,d\eta$ is uniform over the element, and the transformed integral is $$g = \int_{-1}^{1}\!\!\int_{-1}^{1} (4+2\xi)\Big[(4+2\xi)^{2} + (4+3\eta)\Big]\,(6)\,d\xi\,d\eta$$
3. Choose the integration order. The mapped integrand is cubic in $\xi$ and linear in $\eta$. An $n$-point Gauss–Legendre rule integrates polynomials of degree $\le 2n-1$ exactly, so $2n-1 \ge 3$ requires $n = 2$: a $2\times2$ rule with abscissae $\xi_i,\eta_j = \pm 1/\sqrt{3} = \pm 0.5773503$ and unit weights $W_i = W_j = 1$.
4. Locate the stations in physical coordinates. Substituting the abscissae into the map, $$x = 4 \pm 2/\sqrt{3} = 2.845299 \ \text{or}\ 5.154701,\qquad y = 4 \pm 3/\sqrt{3} = 2.267949 \ \text{or}\ 5.732051$$
5. Evaluate in vector–matrix form. Writing the weight vector $\{W\} = \{1\ \ 1\}^{T}$ and collecting the scaled integrand values in the matrix $$[\Phi] = |J|\,\Big[\,f\big(x_i,y_j\big)\,\Big] = \begin{bmatrix} 176.9266 & 236.0650\\ 891.9350 & 999.0734\end{bmatrix}$$ (rows indexed by $\xi_i$, columns by $\eta_j$), the quadrature is the triple product $$g_{2\times2} = \{W\}^{T}[\Phi]\{W\}$$ The individual contributions are tabulated below.

Adding the four contributions, $176.9266 + 236.0650 + 891.9350 + 999.0734$, gives

Part (b) asks what separates this from the exact value, and the honest answer is: nothing. The four-point rule reproduces $g = 2304$ to every digit carried, and two independent conditions conspire to make that happen. First, the map of a rectangle is affine, so $|J|$ is a constant rather than a rational function of $\xi$ and $\eta$; on a general distorted quadrilateral the transformed integrand would not be a polynomial at all and no finite Gauss rule could be exact. Second, the mapped integrand is of degree 3 in $\xi$ and degree 1 in $\eta$, both within the $2n-1 = 3$ exactness limit of the two-point rule in each direction. The quadrature error is therefore identically zero, not merely small.

The instructive contrast is what a cheaper rule delivers. A single-point (centroid) rule samples $f$ at $x = 4$, $y = 4$ with weight $W = 2$ in each direction:

an error of $-16.7\%$. The one-point rule is exact only through degree 1, so it is blind to the curvature of the $x^{3}$ term, and that blindness is exactly the missing sixth. Going the other way, a $3\times3$ rule would also return $2304$, at more than twice the work per element for no gain. The practical lesson the question is testing is that the right quadrature order is the lowest one exact for the integrand at hand — anything lower is wrong, anything higher is wasted. One further check worth carrying: had the Jacobian been omitted, the sum would have come out as $2304/6 = 384$, a factor-of-six error that is easy to spot precisely because $|J|$ is constant.