Evaluating a multiple integral involves expressing it as an iterated integral, which can then be evaluated either symbolically or numerically. We begin by discussing the evaluation of iterated integrals.

firstint=int('x*y','y','1-x','1-x^2')

1/2*x*((1-x^2)^2-(1-x)^2)

1/24

If we begin by declaring x and y to be symbolic variables, we may omit the quotation marks; you will probably want to make use of this most of the time.

firstint=int(x*y,y,1-x,1-x^2)

1/2*x*((1-x^2)^2-(1-x)^2)

1/24

We can even perform the two integrations in a single step:

int(int(x*y,y,1-x,1-x^2),x,0,1)

1/24

There is, of course, no need to evaluate such a simple
integral numerically. However, if we change the integrand to, say, _{}, then MATLAB will be
unable to evaluate the integral symbolically, although it can express the
result of the first integration in terms of erf(x), which is the antiderivative
of _{}.

int(int(exp(x^2-y^2),y,1-x,1-x^2),x,0,1)

Warning: Explicit integral could not be found.

> In D:\matlabr11\toolbox\symbolic\@sym\int.m at line 58

ans =

int(-1/2*erf(-1+x^2)*exp(x^2)*pi^(1/2)+1/2*erf(x-1)*exp(x^2)*pi^(1/2),x = 0 .. 1)

However, we can evaluate the integral numerically:

double(int(int(exp(x^2-y^2),y,1-x,1-x^2),x,0,1))

Warning: Explicit integral could not be found.

> In D:\matlabr11\toolbox\symbolic\@sym\int.m at line 58

ans =

0.1675

** **

** **

We will now address the problem of determining limits for a
double integral from a geometric description of the region of integration.
While MATLAB cannot do that for us, it can provide some guidance through its
graphics and can also confirm that the limits we have chosen define the region
we intended. For a first example, we will consider the bounded region between
the curves y = 2x+2 and y = e^{x}.
We begin by plotting the two curves on the same axes. You may need to
experiment with the interval to get a useful plot; it should be large enough to
show the region of interest, but small enough so that the region of interest
occupies most of the plot.

ezplot(y1,[-1,2]);

hold on;

ezplot(y2,[-1,2]);

hold off; title('region bounded by **y = 2x+2 and y = e^x')** ** **

We can see from the plot that, to define the bounded region
between the two graphs, e^{x}

should be the lower limit for y, and 2x+2 should be the
upper limit. In order to the limits for x are the values for which the two
functions coincide. We use **solve** to find them.

[
-lambertw(-1/2*exp(-1))-1]

[ -lambertw(-1,-1/2*exp(-1))-1]

Hmm…, not clear what this means, but we can convert to numerical values:

-0.7680

1.6783

Now that we have found the limits, we can use **verticalRegion**
to visualize the region they define so that we can check that it is the same as
the region we had in mind. This program
plots the two curves in two different colors: red
for the function entered first, and blue for
the function entered second, and shades the region in between in cyan.

verticalRegion(lims(1),lims(2),y1,y2)

Evidently our limits define the right region. We can now use
them to integrate any function we like over the region in question. To find the
area of the region, for example, we integrate the function 1. Since our limits for *x* are numerical,
a symbolic calculation is not of much use directly, so we use **double** to
convert to a numerical answer.

area=double(int(int(1,y1,y2),lims(1),lims(2)))

2.2270

Of course, the area between two curves could just as easily
have been evaluated as the integral of the difference between the two functions
of *x* whose graphs define the region; that is to say the first
integration in this case is so simple that one can write down the result as
easily as the integral. However, we can also integrate any function we like
over the same region, changing only the integrand. For instance:

newint=double(int(int(sin(x^2+y^2),y,y1,y2),x,lims(1),lims(2)))

Warning: Explicit integral could not be found.

> In C:\MATLABR11\toolbox\symbolic\@sym\int.m at line 58

newint =

0.2892

** **

The numerical evaluation in this case was performed using
numerical routines borrowed from Maple. MATLAB has its own double integrator,
called **dblquad**. There is a complication in using **dblquad**; it does
not accept variable limits. In order to use **dblquad**, we must make a
change of variables in the inner integral. Fortunately, it is always the same
one. We set y=y1+u(y2-y1) so that dy=(y2-y1)du. This gives the identity_{}. **dblquad** demands an inline-vectorized version of the
integrand on the right. All this is easily implemented.

g=inline(vectorize(subs(sin(x^2+y^2),y,y1+u*(y2-y1))*(y2-y1)))

Inline function:

g(u,x) = sin(x.^2+(exp(x)+u.*(2.*x+2-exp(x))).^2).*(2.*x+2-exp(x))

newint2=dblquad(g,0,1,lims(1),lims(2))

Warning: Recursion level limit reached in quad. Singularity likely.

> In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line
102

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In C:\MATLABR11\toolbox\matlab\funfun\quad.m
(quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72

In
C:\MATLABR11\toolbox\matlab\funfun\@inline\feval.m at line 15

In
C:\MATLABR11\toolbox\matlab\funfun\dblquad.m at line 60

Warning: Recursion level limit reached 64 times.

> In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80

In
C:\MATLABR11\toolbox\matlab\funfun\@inline\feval.m at line 15

In
C:\MATLABR11\toolbox\matlab\funfun\dblquad.m at line 60

newint2 =

0.2892

-9.3921e-007

Despite the warnings, the answers provided by these two distinct methods agree to seven decimal places.

Because of the syntactical complications of using **dblquad**
on a region that is not rectangular, we have created **symint2** and **numint2**, which have
identical calling syntax and implement, respectively, symbolic double
integration using **int**, and numerical double integration using **dblquad**.
The advantage of the identical syntax is that if **symint2** fails, one can
try **numint2** instead with a minimum of editing. We illustrate them by
repeating the two previous computations.

double(symint2(sin(x^2+y^2),y,y1,y2,x,lims(1),lims(2)))

Warning: Explicit integral could not be found.

> In C:\MATLABR11\toolbox\symbolic\@sym\int.m at line 58

In
E:\class\math241\common\symint2.m at line 2

ans =

0.2892

numint2(sin(x^2+y^2),y,y1,y2,x,lims(1),lims(2))

Warning: Recursion level limit reached in quad. Singularity likely.

> In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line
102

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In C:\MATLABR11\toolbox\matlab\funfun\quad.m
(quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 142

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141

In
C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72

In
C:\MATLABR11\toolbox\matlab\funfun\@inline\feval.m at line 15

In
C:\MATLABR11\toolbox\matlab\funfun\dblquad.m at line 60

In
E:\class\math241\common\numint2.m at line 7

Warning: Recursion level limit reached 64 times.

> In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80

In
C:\MATLABR11\toolbox\matlab\funfun\@inline\feval.m at line 15

In
C:\MATLABR11\toolbox\matlab\funfun\dblquad.m at line 60

In
E:\class\math241\common\numint2.m at line 7

ans =

0.2892

**Note:** As this example illustrates,
in some problems, **dblquad** and **numint2** do not work that well,
and spit out
long lists of warning messages. In
next week's session, we will discuss
another package called **nit**,
or "numerical integration toolbox", that usually corrects the problem.

a)
Determine the limits of integration for a double integral over
*R*; confirm your limits by using **verticalRegion**.

b)
Evaluate_{}, both using **double** and **int** (or **double**
and **symint2**), and using **dblquad** (or **numint2**).

Let us now turn to the process of integrating over a region
bounded by a level curve of a function of two variables. Once again, the real issue
is the *determination of the limits of integration*. We will consider the region bounded by the
ellipse 3x^{2}+4y^{2}=37. Although it is not strictly necessary
in this case, we will begin by plotting the ellipse, using **genplot**.

genplot(3*x^2+4*y^2,-4:.05:4,-4:.05:4,'contour',[37,37],'k');

axis equal

We can see from the plot that this region is vertically simple. We obtain our y limits by solving the equation of the ellipse for y in terms of x.

[ 1/2*(-3*x^2+37)^(1/2)]

[ -1/2*(-3*x^2+37)^(1/2)]

The two solutions give the upper and lower limits for *y*;
we can find the limits for *x* by setting their difference to 0 and
solving for *x*.

y2=ylims(1)

xlims=solve(y1-y2,x)

-1/2*(-3*x^2+37)^(1/2)

y2 =

1/2*(-3*x^2+37)^(1/2)

xlims =

[ 1/3*111^(1/2)]

[ -1/3*111^(1/2)]

x2=xlims(1)

-1/3*111^(1/2)

x2 =

1/3*111^(1/2)

We can confirm our limits, using **verticalRegion**. The
first two arguments to **verticalRegion** must be numerical, rather than
symbolic. Therefore, we must apply **double** to x1 and x2. verticalRegion(double(x1),double(x2),y1,y2)

We can now integrate any function we desire over the region bounded by the ellipse.

a)
Over the region bounded by the ellipse 3x^{2}+4y^{2}=37.

b)
Over the region bounded by the curve x^{4}+y^{4}=20.

We conclude this lesson with a discussion of integration in
polar coordinates. We consider integrating the function exp(-x^{2}-y^{2})
over the disk bounded by the circle (x-1)^{2 }+ y^{2}=1.
We begin by expressing the function whose vanishing defines the boundary curve
in polar coordinates.

polarfun=simplify(subs((x-1)^2+y^2-1, [x,y],[r*cos(th),r*sin(th)]))

-2*r*cos(th)+r^2

Next we solve the equation of the boundary curve for r in
terms of _{}, to obtain the limits for r.

[ 0]

[ 2*cos(th)]

Since the disk in which we are interested lies entirely to
the right of the x-axis, the appropriate limits for theta are -p/2
and p/2. Bearing in mind that we must include a
factor of *r* when integrating in polar coordinates, we can now set up and
evaluate an iterated integral in any of the usual ways. For example

double(int(int(r*exp(-r^2),0,2*cos(th)),-pi,pi))

Warning: Explicit integral could not be found.

> In C:\MATLABR11\toolbox\symbolic\@sym\int.m at line 58

ans =

2.1724

1)
Evaluate _{}where *R* is the
triangle with vertices at (0,0), (0,2) and (3,0). Hint: begin by determining the equation of the oblique side of
the triangle. Use **verticalRegion** to confirm your limits.

2)
Evaluate _{} where *R* is the region bounded by the closed curve
_{}.

3)
Evaluate _{} where *R *is the region bounded by the cardioid _{}.

4)
Compute the area of the region bounded by the two
curves *x*^{2} - *y*^{2} = 1 and
*y* = *x*^{2} - 3*x* + 1.