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.
We need functions for top and bottom boundaries, according to the integral structure
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.