Skip to main content

Section 18.1 Double integrals over rectangles

A rectangle \(R = [a, b]\times [c, d]\) is the easiest two-dimensional domain to integrate over: in Calculus we would compute

\begin{equation*} \iint_R f(x, y)\,dA = \int_a^b \int_c^d f(x, y)\,dy\,dx \end{equation*}

Suppose we pick two integration rules: one for the interval \([a, b]\text{,}\) say \(\int_a^b g(x)\,dx \approx \sum_i w_i g(x_i)\text{,}\) the other for interval \([c, d]\text{,}\) say \(\int_c^d g(y)\,dy \approx \sum_j v_j g(y_j)\text{.}\) The function \(g\) here is a placeholder; the rule consists of the choice of evaluation points and their weights. Note that the weights and points depend on the interval, so they will probably be different even if we use the same integration method for both variables. With this setup, the two-dimensional integration rule for this rectangle can be expressed as a double sum:

\begin{equation} \iint_R f(x, y)\,dA = \sum_i \sum_j w_i v_j f(x_i, y_j)\label{eq-double-quad-gen}\tag{18.1.1} \end{equation}

Write out explicitly the sum (18.1.1) for the integral of function \(f\) over the rectangle \(R = [-2,2]\times [-1,1]\) using the simple Simpson's rule in both variables.

Answer

In the \(y\) direction, the evaluation points are \(-1, 0, 1\) with the weights \(1/3, 4/3, 1/3\text{.}\) In the \(x\) direction, the evaluation points are \(-2, 0, 2\) with the weights \(2/3, 8/3, 2/3\text{.}\) So the sum has 9 terms:

\begin{equation*} \iint_R f(x, y)\,dA \approx \frac{2}{9} (f(-2,-1) + 4f(0,-1) + f(2,-1) + 4f(-2,0) + 16f(0,0) + 4f(2,0) + f(-2,1) + 4f(0,1) + f(2,1)) \end{equation*}

Given the large number of terms in double sum (18.1.1), we should try to organize such computations efficiently. The first thing to do is to evaluate the function on a rectangular grid formed by all \((x_i, y_j)\) points. A useful Matlab command is meshgrid, which is used as

[X, Y] = meshgrid(x, y)

where x,y are vectors with evaluation points. The output consists of two matrices X, Y which represent the rectangular grid: X has the x-coordinates of all the grid points, and Y has their y-coordinates. The reason this is useful is that f(X, Y) can now be used to evaluate the function on the entire grid at once (producing a matrix of its values), provided \(f\) is written so that it can handle matrix input.

Having computed V = f(X, Y), we need to combine it with weights and to get a single number at the end. Recall in one dimension we would do this as f(x)*w where w is a column vector of x-weights. We should still do this, but to also multiply and sum over the y-weights, we need v*f(X, Y)*w where v is a row of vector of y-weights.

Write a Matlab script to integrate the function \(f(x) = \sqrt{4+xy+y^3}\) over the rectangle in Example 18.1.1, using the integration points and weights from that example.

Answer
x = [-2 0 2];
y = [-1 0 1];
w = (2/3)*[1 4 1];
v = (1/3)*[1 4 1];

f = @(x, y) sqrt(x.*y + y.^3 + 4);
[X, Y] = meshgrid(x, y);
disp(v*f(X, Y)*w')

The first four lines do not involve the function; they are setting the stage. The function is evaluated on the rectangular grid and its values are conveniently aggregated using matrix-vector multiplication. The result is 15.8859. WolframAlpha gives 15.9219.

The rectangular grid computation can also be used to quickly plot the function: surf(X, Y, f(X, Y)) does this.