Evaluating Double Integrals

copyright © 2000 by Paul Green and Jonathan Rosenberg

Evaluating Iterated Integrals

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.

Example 1:

We evaluate the iterated integral . To evaluate the integral symbolically, we proceed in two stages.

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

firstint =

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

answer=int(firstint,'x',0,1)

answer =

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.

syms x y

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

firstint =

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

answer=int(firstint,x,0,1)

answer =

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)

ans =

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

Problem 1

Evaluate the iterated integrals and .

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.

y1=exp(x);y2=2*x+2;

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.

solve(y1-y2)

ans =

[ -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:

lims=double(ans)

lims =

-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)))

area =

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.

syms u

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

g =

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

newint-newint2

ans =

-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.

Problem 2: Let R be the region bounded by the curves y=3x² and y=2x+3.

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).

Integrating over Implicitly Defined Regions

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²+4y²=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.

ylims=solve(3*x^2+4*y^2-37,y)

ylims =

[ 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.

y1=ylims(2)

y2=ylims(1)

xlims=solve(y1-y2,x)

y1 =

-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)]

x1=xlims(2)

x2=xlims(1)

x1 =

-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.

Problem 3:

Integrate

a) Over the region bounded by the ellipse 3x²+4y²=37.

b) Over the region bounded by the curve x⁴+y⁴=20.

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

syms r th

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

polarfun =

-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.

rlims=solve(polarfun,r)

rlims =

[ 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

Problem 4.

Reevaluate the integral above using numint2.

Additional Problems

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² - y² = 1 and y = x² - 3x + 1.