Double integrals over general regions

Section 18.2 Double integrals over general regions

One way to evaluate \(\iint_D f(x, y)\,dA\) for a general region \(D\) is to define \(f\) to be \(0\) outside of \(D\text{,}\) and integrate this redefined function over some rectangle containing \(D\text{.}\) This approach is problematic because the redefined function will probably be discontinuous on the boundary of \(D\text{,}\) and discontinuities cause large errors in numerical integration rules.

🔗

Instead, one can literally adopt the iterated integration approach from calculus courses: do integration in one variable at a time. Suppose the integration is done in \(y\) variable first (“vertical slicing”). Then the limits of the inner integral will be functions of \(x\text{.}\) These functions will need to be implemented in Matlab, and will be called in a loop.

🔗

Example 18.2.1. Implement Simpson’s rule over a non-rectangular region.

Rewrite Example 18.1.2 to integrate the same function over the interior of the ellipse inscribed in the rectangle \(R = [-2,2]\times [-1,1]\text{.}\) (The ellipse with the equation \(x^2/4 + y^2 = 1\text{.}\)) Use composite Simpson’s rule with \(n=4\) subintervals in both variables.

🔗

Answer.

We need functions for top and bottom boundaries, according to the integral structure

\begin{equation*} \int_a^b \int_{\text{bottom}}^{\text{top}} f(x, y)\,dy\,dx \end{equation*}

In the x-variable, the points and weights are set once. But in the y-direction they are calculated within a loop, because different x-values produce different limits of integration for y. The vector innerint is used below to store the values of the inner integral.

🔗

f = @(x, y) sqrt(x.*y + y.^3 + 4);
top = @(x) sqrt(1-x.^2/4);
bottom = @(x) -sqrt(1-x.^2/4);
a = -2;
b = 2;

n = 4;
x = linspace(a, b, n+1);
w = [1 4 2 4 1]*(b-a)/(3*n);
innerint = zeros(size(x));

for k = 1:numel(x)
    c = bottom(x(k));
    d = top(x(k));
    y = linspace(c, d, n+1)';
    v = [1 4 2 4 1]*(d-c)/(3*n);
    innerint(k) = v*f(x(k), y);
end
    
disp(innerint*w');

The result is 11.8764, compared to the WolframAlpha output 12.5419. This is less accurate than what we got for rectangle. Note that all points with \(x=\pm 2\) are wasted because there is no interval to integrate over with respect to \(y\text{.}\) For this reason, open-ended rules such as Gaussian integration is preferable. Replacing the choice of x, w above with the 3-point Gaussian rule

x = [-sqrt(3/5) 0 sqrt(3/5)]*(b-a)/2 + (a+b)/2;
w = [5/9 8/9 5/9]*(b-a)/2;

we get 12.7081, which is much more accurate.

🔗