Processing math: 100%
Skip to main content
🔗

Section 18.2 Double integrals over general regions

🔗
One way to evaluate ∬Df(x,y)dA for a general region D is to define f to be 0 outside of D, and integrate this redefined function over some rectangle containing D. This approach is problematic because the redefined function will probably be discontinuous on the boundary of D, 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. These functions will need to be implemented in Matlab, and will be called in a loop.

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.