Chapter 12	Vector Geometry

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

12.1  Vectors

Students should read Sections 12.1 - 12.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

A vector is an object that has magnitude and direction.  In physics, these vectors are denoted by arrows, where the magnitude of the vector is represented by the length of the vector, and the way in which the arrow points indicates its direction.  In mathematics, these vectors are represented by points in two or three dimensions, where the vector is the arrow that starts at the origin and ends at the point.  For example, the point (2, 1, 3) could be considered both as a point in 3-D space and as a vector from (0, 0, 0) to (2, 1, 3). To distinguish a point from a vector, we will use the angled brackets ⟨ and ⟩ instead of parentheses. Thus, the point (2, 1, 3) is denoted (2, 1, 3) as usual, but the vector from the origin to that point is denoted ⟨2, 1, 3⟩.

The length or magnitude of a vector v is denoted ||v||, and is read as "norm v."  If  v = ⟨a, b, c⟩, then ||v|| = (Sqrt[(a^2) + (b^2) + (c^2)]).  In two dimensions, if v = ⟨a, b⟩, then ||v|| = (Sqrt[(a^2) + (b^2)]).

Vectors and matrices, in Mathematica, are simply lists.  A vector is a list of numbers within braces, with commas between numbers, while a matrix is a list of lists (vectors), with each vector list being a row of the matrix (for a complete description of lists in Mathematica, see Section 1.2.3 of this text).  Of course, all rows must be the same size.  For example, consider the vector a below:

a = {1,3 ,5}

The ith component of the vector a is denoted by (Subscript[a, i]), or in Mathematica, by a[[i]].  For instance the second component of a, which is 3, would be obtained by:

a[[2]]

All of the usual vector algebra operations are available to us:

Dot Product

The Dot Product of two vectors u=⟨(Subscript[u, 1,]) (Subscript[u, 2]), (Subscript[u, 3])⟩ and v=⟨(Subscript[v, 1,]) (Subscript[v, 2]), (Subscript[v, 3])⟩ is defined by 
u·v=(Subscript[u, 1])(Subscript[v, 1]) + (Subscript[u, 2])(Subscript[v, 2]) + (Subscript[u, 3])(Subscript[v, 3]).  For example:

a = {1, 3, 5}
b = {1, -2, 3}
a.b

or

Dot[a,b]

NOTE: We use the ordinary period symbol on the keyboard for the dot product.

Cross Product

The cross product of two vectors u=⟨(Subscript[u, 1,]) (Subscript[u, 2]), (Subscript[u, 3])⟩ and v=⟨(Subscript[v, 1,]) (Subscript[v, 2]), (Subscript[v, 3])⟩, is defined as a vector perpendicular to both u and v, and calculated by the following "right-hand" rule:

u × v = ⟨(Subscript[u, 2])(Subscript[v, 3]) - (Subscript[u, 3])(Subscript[v, 2]), (Subscript[u, 3])(Subscript[v, 1]) -  (Subscript[u, 1])(Subscript[v, 3]) , (Subscript[u, 1])(Subscript[v, 2]) - (Subscript[u, 2])(Subscript[v, 1])⟩

This calculation can be done in Mathematica in two ways.  The first is to use the Cross command:

Cross[a,b]

The second is by using the multiplication symbol "×".  This special symbol can be entered on the Basic Math Input Palette or by pushing the escape key, followed by typing the word "cross" and hitting the escape key again:  [esc]cross[esc]

a \[Cross] b

Recall that the cross product of 2 vectors, a and b creates a vector perpendicular to the plane of the vectors a and b.  In your Calculus text, the cross product is also defined as the determinant of a special matrix. We will look at this a little later.

Norm (Length) of a Vector

The norm or length of a vector can be calculated in Mathematica by the Norm command

Clear[x,y,z]

Norm[{x,y,z}]

Norm[a]

Norm [2a]

Vector Addition

The sum of two vectors u=⟨(Subscript[u, 1,]) (Subscript[u, 2]), (Subscript[u, 3])⟩ and v=⟨(Subscript[v, 1,]) (Subscript[v, 2]), (Subscript[v, 3])⟩ is defined to be u+v=(Subscript[u, 1])(Subscript[v, 1]) + (Subscript[u, 2])(Subscript[v, 2]) + (Subscript[u, 3])(Subscript[v, 3]).

2a - 3b +{1,1,1}

Example 12.1.  Let a=⟨1,2,3⟩.  Show that (a/(||a||))  is a unit vector.

Solution:

Norm[a/Norm[a]]

Example 12.2.  Find the equation of a line in 3-space passing through (Subscript[P, 0]) = (3,-1,4) in the direction of  v = ⟨2,7,1⟩ and graph it.

Solution: The line through (Subscript[P, 0]) =((Subscript[x, 0]),(Subscript[y, 0,])(Subscript[z, 0])) in the direction of v=⟨a,b,c⟩ is described in vector or parametric form by:

Vector form:  r(t)=⟨(Subscript[x, 0]),(Subscript[y, 0,])(Subscript[z, 0])⟩+t⟨a,b,c⟩
Parametric Form:  x=(Subscript[x, 0])+a t, y=(Subscript[y, 0])+b t, z=(Subscript[z, 0])+c t

Thus, the vector description of the line is

Clear[r,t];
r[t_]= {3,-1,4} + t {2,7,1}

To graph this line we use the ParametricPlot3D command:

ParametricPlot3D[{(Subscript[f, x]),(Subscript[f, y]),(Subscript[f, z])},{u,(Subscript[u, min]),(Subscript[u, max])}]
produces a three-dimensional space curve parametrized by a variable u which runs from (Subscript[u, min]) to (Subscript[u, max]).

ParametricPlot3D[r[t], {t,-3,3},ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

NOTE: This plot command uses the option ImageSize to specify the size of graphics output.  Settings include Tiny, Small, Medium, Large, or {pt}, where pt is the number of points.

Example 12.3.  Give the description in vector form of the line that passes through the points P=(1,0,4) and Q=(3,2,1), then find the midpoint of the line segment (OverVector[PQ]) and plot this line segment.

Solution: The line through points P=((Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1])) and Q=((Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2])) has vector form r(t)=(1-t) ⟨(Subscript[a, 1]),(Subscript[b, 1]),(Subscript[c, 1])⟩+t ⟨(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2])⟩.  In this parametrization, r(0) = P and r(1) = Q.  Thus,

r[t_] = (1-t) {1,0,4} + t {3,2,1}

The midpoint of the line segment (OverVector[PQ]) is

r[(1/2)]

The plot of the line segment is

ParametricPlot3D[r[t], {t,-0.1,1.1},ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

Example 12.4.  Find the angle between the vectors v = ⟨3,6,2⟩ and w = ⟨6,3,6⟩.

Solution: Remember that the angle between two vectors, v and w, is given by θ, which is defined by θ=(cos^(-1))(((v.w)/(||v|| ||w||))) .  Therefore,

v= {3,6,2}
w = {6,3,6}

θ = ArcCos[((v.w)/(Norm[v]Norm[w]))]

N[%]

Therefore, θ = .7045 radians.

Exercises

1. Calculate the length of the vector v=⟨1,3,4⟩.

In Exerices 2 and 3, calculate the linear combinations.
2.   5⟨2,-2,5⟩+6⟨1,3,8⟩ 		3.   6 ⟨2, 0, -1⟩  - 3 ⟨8, 6, 9⟩

4. Find a vector parametrization for the line that passes through P=(1,2,-6) with direction vector v= ⟨2,1,5⟩.

In Exercises 5 and 6, determine whether the two given vectors are orthogonal (v ⊥ w iff  v.w = 0):
5.  ⟨1,1,1⟩,  ⟨1,-2, 3⟩		6.  ⟨1,1,1⟩,  ⟨-3,2,1⟩

In Exercises 7 and 8, find the angle between the vectors:
7.  ⟨1,2⟩,  ⟨5,7⟩			8.  ⟨2,4,1⟩,  ⟨1,-3,5⟩

12.2  Matrices and the Cross Product

Students should read Section 12.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

In order to understand the alternate approach to the cross product alluded to above, we need to define the terms matrix and determinant.

Matrices

A matrix is a rectangular array of numbers consisting of n rows and m columns (denoted n × m).  We are especially interested in square matrices where m = n and, in particular, m = 2 or m=3.  For example:  A 3 × 3 matrix would be

(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), (Subscript[a, 33])}}))

but Mathematica would show this matrix as:

A = Table[10i+j,{i,3},{j,3}]

B = Table[i+j,{i,2},{j,2}]

To have Mathematica display a matrix in the traditional way, use the MatrixForm command:

MatrixForm[A]
MatrixForm[B]

Note that in the definition of the matrices A and B, Mathematica treats them as lists and when we use the command MatrixForm, we can see the matrices presented in the traditional way.

Determinants

The determinant is a function, Det, which assigns to each square matrix a number which is defined for 2 × 2 and 3 × 3 matrices as follows:

Clear[a,b];
F = {{a,b},{c,d}}
MatrixForm[F]

Det[F]

G= {{a1,a2,a3},{b1,b2,b3},{c1,c2,c3}}
MatrixForm[G]

Det[G]

Using these definitions, we can now define the cross product of two vectors by the formula

⟨(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])⟩×⟨(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3])⟩=Det(({{i, j, k}, {(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])}, {(Subscript[c, 1]), (Subscript[c, 2]), (Subscript[c, 3])}}))

where i =⟨1,0,0⟩, j =⟨0,1,0⟩, and k =⟨0,0,1⟩.

Example 12.5.  Calculate the cross product of  v=⟨1,3,6⟩ and w=⟨-2,8,5⟩.

Solution:

Clear[i,j,k]
g = {i,j,k}
v= {1,3,6}
w = {-2,8,5}
A = {g,v,w}

MatrixForm[A]

v \[Cross] w
Det[A]

Observe that the two previous outputs are equivalent.

Exercises

1.  Calculate the determinants of (({{0, 5, 0}, {1, 3, 6}, {2, 5, 5}})) and of (({{3, 5}, {6, 2}})).

2.  Calculate the cross product of  v=⟨2,0,0⟩ and w=⟨-1,0,1⟩.  Do this using the Cross command as well as by the determinant approach.

3.  Calculate the area of the parallelogram spanned by the vectors v and w above. (Hint: look up the formula for this in your calculus textbook.)

4.  Calculate the volumn of the parallelepiped spanned by:
      u  =  ⟨2, 2, 1⟩,  v  =  ⟨1, 0, 3⟩,   and  w  =  ⟨0, -4, 2⟩

5.  Show that  v×w  = -w×v  and that  v×v = 0.

12.3  Planes in 3-Space

Students should read Section 12.5 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Note that a plane in 3-D space is defined as all points P(x,y,z) such that the line segment (OverVector[(Subscript[P, 0])P ])  is perpendicular to a given vector n, called the normal vector, where the initial point of n is (Subscript[P, 0])=((Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])).  In vector notation, this is described by the equation n· (OverVector[(Subscript[P, 0])P ])  = 0, where  (OverVector[(Subscript[P, 0])P ]) = ⟨x-(Subscript[x, 0,]) y-(Subscript[y, 0]),z-(Subscript[z, 0])⟩.  Therefore, the equation of the plane through (Subscript[P, 0]) = ((Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])) with nonzero normal vector n = ⟨a,b,c⟩ can be denoted by either of the following:

Vector form:			n·⟨x,y,z⟩=d
Scalor form: 			a x + b y + c z = d

Here, d=a (Subscript[x, 0])+b (Subscript[y, 0])+c (Subscript[z, 0])=n·⟨(Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])⟩.

Example 12.6.  Find an equation of the plane determined by the points P=(1,0,-1), Q=(2,2,1), and R=(4,2,5).  Then plot the graph of the plane.

Solution: The vectors a=(OverVector[PQ ]) and  b=(OverVector[PR])  lie in the plane, so the cross product n=a×b is normal to the plane:

Clear[a, b, n]
a = {2, 2, 1} - {1, 0, -1}
b = {4, 2, 5} - {1, 0, -1}
n = a \[Cross] b
n . {x,y,z}==d

To compute the value of d, we choose any point on the plane, that is, we can choose either P, Q, or R, and then compute d=n·P, d=n·Q, or d=n·R.  Let us choose P=⟨1,0,-1⟩.

d = n . {1, 0, -1}

Therefore, the plane we want has equation 8x-4z=12  and the graph is obtained by using the ContourPlot3D command which has the form:

({{, ContourPlot3D[f,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])}]
which produces a three-dimensional contour plot of f as a function of x, y and z. 
or}, {, ContourPlot3D[f==g,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},{z,(Subscript[z, min]),(Subscript[z, max])}]
which plots the contour surface for which f=g. }})

ContourPlot3D[8x - 4z==12,{x,-2,2},{y,-2,2},{z,-2,2},ImageSize->{250}]

In order to see this plane more clearly from a different perspective, move your cursor over the plot.  Then drag the mouse while pressing and holding the left mouse button to rotate the plot.

Exercises

1.  Let PL be the plane with equation 7x-4y+2z=10.  Find an equation of the plane QL parallel to                       PL and passing through Q=(2,1,3) and graph it.

2.  Find the equation of the plane through the points P=(1,5,5), Q=(0,1,1), and R=(2,0,1) and     graph it.

3.  Find the angle between the two planes:  x+2y+z=3 and 4x+y+3z=2.  (Hint: The angle between two planes is the angle between their normal vectors.)

12.4  A Survey of Quadric Surfaces

Students should read Section 12.6 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

A quadric surface is the three-dimensional equivalent of a conic section (i.e., ellipses, hyperbolas, and parabolas).  The basic types of quadric surfaces are ellipsoids, hyperboloids (of one or two sheets), paraboloids (elliptic or hyperbolic), and cones.

12.4.1  Ellipsoids

The standard ellipsoid is described by (((x/a))^2)+(((y/b))^2)+(((z/c))^2)=1. To help us visualize it, we are often interested in the mesh of curves called traces, obtained by intersecting our quadric surface with planes parallel to one of the coordinate planes. In the plot below, you can see that mesh, and also see that the traces of an ellipsoid are themselves ellipses.

Example 12.7.  Graph the ellipsoid above, with a=3, b= 4, and c = 5, and describe the traces of this ellipsoid.

Solution: The correct Mathematica command to use is ContourPlot3D.  This is shown following:

ContourPlot3D[(x/3)^2+(y/4)^2+(z/5)^2==1,{x,-6,6},{y,-6,6},{z,-6,6},AxesLabel->{x,y,z},ImageSize->{250}]

Again, note that the ellipsoid can be manually rotated to look at it from different perspectives. First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.  When you do this, you will note that, indeed, all of the traces are ellipses.

12.4.2  Hyperboloids

The three-dimensional hyperbolas are called hyperboloids, and come in two types: the hyperboloid of one sheet, with standard form (((x/a))^2)+(((y/b))^2)=(((z/c))^2)+1, and the hyperboloid of two sheets, with standard form (((x/a))^2)+(((y/b))^2)=(((z/c))^2)-1.  A limiting case of the hyperboloid is the elliptic cone, defined by the equation (((x/a))^2)+(((y/b))^2)=(((z/c))^2).

Example 12.8.  Describe the traces of the two hyperboloids: (((x/3))^2)+(((y/4))^2)=(((z/5))^2)+1 and (((x/3))^2)+(((y/4))^2)=(((z/5))^2)-1.

Solution: First we graph the hyperboloids:

ContourPlot3D[(x/3)^2+(y/4)^2   == (z/5)^2+ 1,{x,-6,6},{y,-6,6},{z,-6,6},AxesLabel->{x,y,z},ImageSize->{250}]

In this case, the traces parallel with the xy-axis are all ellipses, and the traces parallel wth the xz- and yz-axes are hyperbolas.

ContourPlot3D[(x/3)^2+(y/4)^2 ==(z/5)^2-1,{x,-30,30},{y,-30,30},{z,-30,30},AxesLabel->{x,y,z},ImageSize->{250}]

When we look at this plot, we see that the traces are the same as for the previous hyperboloid of one sheet.

Example 12.9.  Graph the cone with a=3, b=4, and c=5, and define its relationship to the hyperboloid of one sheet.

Solution: We get the graph by using the ContourPlot3D Command:

ContourPlot3D[(x/3)^2+(y/4)^2 ==(z/5)^2,{x,-30,30},{y,-30,30},{z,-30,30},AxesLabel->{x,y,z},ImageSize->{250}]

When we compare this plot with that of the hyperboloid of one sheet (see previous example), we can see clearly that this cone can be thought of as a limiting case of the hyperboloid of one sheet in which we pinch the waist down to a point.

12.4.3  Paraboloids

The final family of quadric surfaces that we want to consider are the paraboloids, of which there are two types: elliptic and hyperbolic. Their standard equations are z=(((x/a))^2)+(((y/b))^2) (elliptic paraboloid) and z=(((x/a))^2)-(((y/b))^2) (hyperbolic paraboloid).

Example 12.10.  Graph the two types of paraboloids for a=3 and b=4 and describe their traces.

Solution: Here is the graph of the elliptic paraboloid:

ContourPlot3D[(x/3)^2+(y/4)^2 ==z,{x,-30,30},{y,-30,30},{z,-30,30},AxesLabel->{x,y,z},ImageSize->{250}]

Observe that the traces in the direction of the xz- and yz-axes are both parabolas while those in the xy-direction are ellipses, which can be seen by dragging the plot in the appropriate directions.  Similarly, for the hyperbolic paraboloid:

ContourPlot3D[(x/3)^2-(y/4)^2 ==z,{x,-30,30},{y,-30,30},{z,-30,30},AxesLabel->{x,y,z},ImageSize->{250}]

Again, by dragging the plot above, we see that the traces in the yz-direction are parabolas while those in the xz-direction are hyperbolas.

12.4.4  Quadratic Cylinders

The last group of quadric surfaces we will look at are the quadratic cylinders. These are surfaces formed from a two-dimensional curve (in the xy-plane) along with all vertical lines passing through the curve:

Example 12.11.  Graph a selection of quadratic cylinders.

Solution:

a)  A circular cylinder of radius r:  (x^2)+(y^2)=(r^2).  For the graph, we will use r=3.

ContourPlot3D[x^2+y^2 ==3^2,{x,-5,5},{y,-5,5},{z,-30,30},AxesLabel->{x,y,z},ImageSize->{250}]

b)  An elliptic with equation (((x/a))^2)+(((y/b))^2)=1. We will use a=3 and b=6.

ContourPlot3D[(x/3)^2+(y/6)^2 ==1,{x,-5,5},{y,-8,8},{z,-20,20},AxesLabel->{x,y,z},ImageSize->{250}]

c) A hyperbolic cylinder with equation (((x/a))^2)-(((y/b))^2)=1.  We will use a=3 and b=6.

ContourPlot3D[(x/3)^2-(y/6)^2 ==1,{x,-10,10},{y,-10,10},{z,-20,20},AxesLabel->{x,y,z},ImageSize->{250}]

d) A parabolic cylinder with equation y =a (x^2) with a=3.

ContourPlot3D[y ==3x^2,{x,-3,3},{y,-1,8},{z,-10,10},AxesLabel->{x,y,z},ImageSize->{250}]

Exercises

In Exercises 1 through 5, state whether the given equation defines an ellipsoid, hyperboloid, or paraboloid, and of which type.  Then confirm your answer by graphing the quadric surface.
1.   (((x/5))^2)+(((y/7))^2)+(((z/9))^2)=1
2.   (((x/5))^2)-(((y/7))^2)+(((z/9))^2)=1
3.    (x^2)+5(y^2)-6(z^2)=1
4.    z=(((x/5))^2)+(((y/7))^2)
5.    z=(((x/5))^2)-(((y/7))^2)

In Exercises 6 through 8, state the type of  the quadric surface and graph it, and then describe the trace obtained by intersecting it with the given plane.
6.   (((x/5))^2)+(y^2)+(((z/9))^2)=1,  z=1/4
7.   y=2(x^2),    z=25
8.   (((x/5))^2)-(((y/7))^2)+(((z/9))^2)=1,  y=4

12.5  Cylindrical and Spherical Coordinates

Students should read Section 12.7 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

12.5.1  Cylindrical Coordinates

In cylindrical coordinates, the point P=(x,y,z) is expressed as (r,θ,z) where r and θ are the polar coordinates of x and y.  The formulas for converting from (x,y,z) to (r,θ,z) are:

Cylindrical to Rectangular		Rectangular to Cylindrical
x=r cos θ				r=(Sqrt[(x^2)+ (y^2)])
y=r sin θ				tan θ=y/x
z=z					z=z

The commands in Mathematica to do these conversions must first be loaded into Mathematica from the "Vector Analysis" external package:

<<VectorAnalysis`

Example 12.12.  Convert (r,θ,z)=(2,3π/4,5) to rectangular coordinates.

Solution: We use the CoordinatesToCartesian command to convert from cylindrical to rectangular coordinates:

CoordinatesToCartesian[{2,3Pi/4,5}, Cylindrical]

N[%]

Example 12.13.  Convert (x,y,z)=(2,3,5) to cyclindrical coordinates.

Solution: We use the CoordinatesFromCartesian command to convert from rectangular to cylindrical coordinates:

CoordinatesFromCartesian[{2, 3, 5},Cylindrical]

N[%]

Of course, one very strong point for Mathematica is its graphing ability.  It will easily graph functions described in cylindrical coordinates. The command to do this is RevolutionPlot3D.

RevolutionPlot3D[(Subscript[f, z]),{t,(Subscript[t, min]),(Subscript[t, max])},{θ,(Subscript[θ, min]),(Subscript[θ, max])}]
takes the azimuthal angle θ to vary between (Subscript[θ, min]) and (Subscript[θ, max]).

Example 12.14.  Graph the cylindrical coordinate function z=((2(r^2)sin(5θ))/(1+(r^2))).

Solution:

Clear[r,θ];
RevolutionPlot3D[((2(r^2) Sin[5 θ])/(1+(r^2))),{r,0,5},{θ,0,2 π},ImageSize->{250}]

12.5.2  Spherical Coordinates

A point P=(x,y,z) is described in spherical coordinates by a triple (ρ,θ,ϕ) where ρ is the distance of P from the origin, θ is the polar angle of the projection (x,y,0), and ϕ is the angle between the z-axis and the ray from the origin through P. The formulas for converting between rectangular and spherical coordinates are:

Spherical to Rectangular		Rectangular to Spherical

x=ρ cos θ sin ϕ			ρ=(Sqrt[(x^2) + (y^2) + (z^2)])
y=ρ sin θ sin ϕ			tan θ=y/x
z=ρ cos ϕ				cos ϕ=z/ρ

These conversions are done in Mathematica using the same commands as with cylindrical coordinates, but with the word spherical replacing cylindrical.

Example 12.15.  Convert (ρ,θ,ϕ)=(2,3π/4,π/5) to rectangular coordinates.

Solution:

CoordinatesToCartesian[{2,3Pi/4,π/5}, Spherical]

N[%]

Example 12.16.  Convert (x,y,z)=(2,3,5) to spherical coordinates.

Solution:

CoordinatesFromCartesian[{2, 3, 5},Spherical]

N[%]

Again, the main use here of Mathematica is its graphing ability.  It will easily graph functions described in spherical coordinates. The command to do this is the SphericalPlot3D command.

SphericalPlot3D[r,{θ,(Subscript[θ, min]),(Subscript[θ, max])},{ϕ,(Subscript[ϕ, min]),(Subscript[ϕ, max])}] 
generates a 3D plot with a spherical radius r as a function of spherical coordinates θ and ϕ.

Example 12.17.  Graph the spherical coordinate function ρ=1+sin (6ϕ)/6.

Solution:

SphericalPlot3D[ ρ = 1+Sin[6ϕ]/6,{θ,0,Pi},{ϕ,0,2Pi},ImageSize->{250}]

Exercises

Convert from cylindrical to rectangular:
1.  (2,π/3,-4)  		2.  (1,π/2,3)

Convert from rectangular to cylindrical:
3.   (2,2,5)			4.   (4,(Sqrt[3]),8)

5.  Plot the surface (z^2)+(r^2)=25θ and describe it.

Convert from spherical to rectangular:
6.    (2,π/5,π/3)		7.    (4,π/6,5π/6)

Convert from rectangular to spherical:
8.   ((Sqrt[2]),2,3)			9.   (4,(Sqrt[3])/2,(Sqrt[8]))

10.  Plot the surface ρ sin ϕ=5 and describe it.

Chapter 13	Calculus of Vector-Valued Functions

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

13.1.  Vector-Valued Functions

Students should read Section 13.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

A vector-valued function is a vector where the components of the vector are themselves functions of a common parameter (or variable). For example, r is a vector-valued function if r(t) = ⟨x(t),y(t),z(t)⟩.  If we think of t as the time variable, the r(t) describes the motion of a particle through three-dimensional space over time.  What we want to do is to understand what path is taken.  We do this through graphing in three dimensions. Also, sometimes it is helpful to consider the projections of these curves onto the coordinate planes.  For example, the projection of r(t) on the xy-plane is ⟨x(t),y(t),0⟩.

Example 13.1. Trace the paths of each of the following vector functions and describe its projections onto the xy-, xz-, and yz-planes:
a)  r(t)=⟨t, (t^2), 2t⟩
b)  r(t)=⟨(cos^3)t, (sin^3)t,sin 2t⟩

Solution:  We use the ParametricPlot3D command to trace the path of each curve and to see its projection.

a) First, we look at the plot of r(t)=⟨t,(t^2), 2t⟩:

ParametricPlot3D[{t,(t^2), 2t},{t,-3,3},PlotStyle-> Red,ImageSize->{250}]

This curve looks very much like a parabola in 3-D space.  To see the projections, we look first at:

ParametricPlot3D[{t,(t^2), 0},{t,-3,3},PlotRange->{-1,1},PlotStyle-> Orange,ImageSize->{250}]

This is clearly a parabola in the xy-plane.

ParametricPlot3D[{t,0, 2t},{t,-3,3},Ticks->{Automatic,{-1,0,1},Automatic},PlotStyle-> Orange,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

And this clearly a line in the xz-plane.

ParametricPlot3D[{0,(t^2), 2t},{t,-3,3},Ticks->{{-1,0,1},Automatic,Automatic},PlotStyle-> Orange,ImageSize->{250}]

This last plot is also clearly a parabola, but in the yz-plane.

b)  Next, we look at r(t)=⟨(cos^3)t ,(sin^3) t,sin 2t⟩:

ParametricPlot3D[{((Cos[t])^3),((Sin[ t])^3),Sin[2t]},{t,-2π,2π},PlotStyle-> Orange,ImageSize->{250}]

Note that since both sine and cosine are periodic with period 2π, it is not necessary to extend the domain beyond -2π or +2π. The projection in the xy-plane is:

ParametricPlot3D[{((Cos[t])^3),((Sin[ t])^3),0},{t,-2π,2π},PlotPoints->100,ImageSize->{250}]

The projection in the xz-plane is:

ParametricPlot3D[{((Cos[t])^3),0,Sin[2t]},{t,-2π,2π},ImageSize->{250}]

Lastly, the projection in the yz-plane is:

ParametricPlot3D[{0,((Sin[ t])^3),Sin[2t]},{t,-2π,2π},ImageSize->{250}]

Note that the last two projections are almost exactly alike.  This is to be expected because the sine and cosine functions have the same graph, but π/2 radians apart.

Exercises

In Exercises 1 through 3,  graph r(t) and its three projections onto the coordinate planes.
1.  r(t)=⟨cos( 2t),cos t,sin t⟩		2.  r(t)=⟨t+15,(e^(0.08t))cos t,(e^(0.08t))sin t⟩
3.  r(t)=⟨t,t,25t/(1+(t^2))⟩

4. Which of the following curves have the same projection onto the xz-plane?  Graph the three projections to check your answer.
a.  (Subscript[r, 1])(t)=⟨t,(e^t),(t^2)⟩		b.  (Subscript[r, 2])(t)=⟨(e^t),t,(t^2)⟩		c.  (Subscript[r, 3])(t)=⟨t,cos t,(t^2)⟩

13.2.  Calculus of Vector-Valued Functions

Students should read Section 13.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Since vector-valued functions are differentiated and integrated component by component, Mathematica will handle this easily since it treats vectors as lists and automatically performs the indicated operation on each element of the list.

The derivative of a vector valued function r(t)=⟨x(t),y(t),z(t)⟩ is defined to be

r'(t)=⟨x'(t),y'(t),z'(t)⟩

while the integral of r(t) is

∫r(t)dt=⟨∫x(t)dt,∫y(t)dt,∫z(t)dt⟩.

Similarly, the limit is defined by

(lim_(t->a))r(t)=⟨(lim_(t->a))x(t),(lim_(t->a))y(t),(lim_(t->a))z(t)⟩.

Example 13.2.  Differentiate and integrate each of the following vector functions:
a) r(t)=⟨t, (t^2), 2t⟩
b) s(t) = ⟨(cos^3) t , (sin^3) t, sin 2t⟩

Solution:

(a)

Clear[r,s,t]

r[t_]:= {t,(t^2),2t}
s[t_]:= {((Cos[t])^3),((Sin[t])^3),Sin[2t]}

(∂_t) r[t]

∫r[t]\[DifferentialD]t

(b)

(∂_t) s[t]

∫s[t]\[DifferentialD]t

Limits are handled the same way both in the calculus of vector-valued functions and in Mathematica:

Example 13.3.  Evaluate (Underscript[limit, h->0]) ((r(t+h)-r(t))/h)  for r(t)=⟨t, (t^2), 2t⟩.

Solution:

Since r(t) has been defined in the previous example, we merely evaluate

Limit[((r[t+h]-r[t])/h),h->0]

As we would expect, this limit gives us the same answer for r'(t) as in the previous example.

Example 13.4.  Evaluate (Underscript[limit, t->3])⟨(t^2),4t,(1/(t^3))⟩.

Solution:

Limit[{(t^2), 4t,(1/(t^3))},t->2]

Derivatives of Dot and Cross Products

Using the formulas of the derivative of the dot and cross products for vector-valued functions is simple in Mathematica. As a reminder, the formulas are:

(d/dt)(r(t)·s(t))=r(t)·s'(t)+r'(t)·s(t) and (d/dt)(r(t)×s(t))=r(t)×s'(t)+r'(t)×s(t)

Example 13.5.  Evaluate (d/dt)(r(t)·s(t)) and (d/dt)(r(t)×s(t)) for r(t)=⟨t,(t^2),2t⟩ and s(t) = ⟨(cos^3) t ,(sin^3) t,sin 2t⟩.

Solution:

(∂_t)(r[t].s[t])

(∂_t)(r[t]\[Cross]s[t])

Tangent Lines

Example 13.6.  Find the vector parametrization of the tangent line to r(t)=⟨1-(t^2),5t,(t^3)⟩ at the point t=1 and plot it along with r(t).

Solution: Recall that the tangent line at (Subscript[t, 0]) has vector parametrization L(t)=r((Subscript[t, 0]))+t r'((Subscript[t, 0])):

r[t_]= {1-(t^2),5t,(t^3)}
r'[t]
L[t_]=r[1]+t*r'[1]

Here is a plot of the curve and the tangent line.

ParametricPlot3D[{r[t],L[t]},{t,-2,4},ImageSize->Small]

NOTE: Recall that the plot can be rotated to better view it from different perspectives.

Exercises

In Exercises 1 and 2 evaluate the limits
1.  (lim_(t->π))⟨sin 2t,cos t,tan 4t⟩ 		2.  (lim_(t->0))⟨(1/(t+1)),(((e^t)-1)/t),4t⟩

In Exercises 3 and 4 compute the derivative and integral.
3.  r(t)=⟨tan t,4t-2,sin t⟩			4.  r(t)=⟨(e^t),(e^(2t))⟩

5.  Find a parametrization of the tangent line at the point indicated and plot both the vector-valued curve and the tangent line on the same set of axes.

6.  Evaluate (d/dt)r(g(t)) for r(t)=⟨4 sin 2t,2 cos 2t⟩ and g(t)=(t^2).

13.3.  Arc Length

Students should read Section 13.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

The arc length of a path r(t)=⟨x(t),y(t),z(t)⟩ for a≤t≤b is given by

L=(∫_a^b)||r'(t)||dt=(∫_a^b)(Sqrt[(((x'(t)))^2) + (((y'(t)))^2) + (((z'(t)))^2)])dt

and like the one-dimensional version is difficult to evaluate by hand.  Thus Mathematica is the perfect tool for calculating this.

Example 13.7.  Compute the arc length of r(t)=⟨1-(t^2),5t,2 (t^3)⟩  over the interval 1≤t≤2.

Solution:

r[t_] := {1-(t^2), 5t,2 (t^3)} 
L = (∫_1^2)Norm[r'[t]] \[DifferentialD]t

Note that the above output indicates that Mathematica cannot find an antiderivative for the integrand, and thus we need to find another technique to evaluate this integral.  Hence, we next try the numerical integrate command, NIntegrate, which does give us our result:

L = NIntegrate[Norm[r'[t]],{t,1,2}]

Speed

The vector r'(t) is also known as the velocity vector as it points in the (instantaneous) direction of motion described by r(t).  Its length or norm, ‖r'(t)‖, gives the speed at time t.

Example 13.8.  Compute the speed of r(t)=⟨1-(t^2),5t,2 (t^3)⟩ when t=1, 1.5, and 2.

Solution:

The following output gives a list of speeds of r'(t) at the three given times using the Norm command, which calculates the norm of a vector:

r[t_] := {1-(t^2), 5t,2 (t^3)} 
Speed = {Norm[r'[1]], Norm[r'[1.5]], Norm[r'[2]]}

N[%]

Observe that the speed is increasing as we move along the path of r(t) from t=1 to t=2.  This can be seen graphically by plotting the speed:

Norm[r'[t]]
Plot[Norm[r'[t]],{t,1,2}]

NOTE: Observe how the Norm command inserts absolute values around each vector component in the formula for ‖r'(t)‖, which seems redundant since each component is squared.  This is done because in Mathematica vector components are allowed to be complex-valued, in which case absolute values are needed to refer to their magnitudes.

Exercises

In Exercises 1 and 2, compute the length of curve over the given interval.
1.  r(t)=⟨2 sin t,6t,2cos t⟩,   -6≤t≤6		2.  r(t)=⟨12 t,8(t^(3/2)),3 (t^2)⟩,   0≤t≤1

In Exercises 3 and 4, find the speed of a particle moving along the curve r(t) at the given value of t.
3.  r(t)=⟨(e^(t-2)),15 t,5/t⟩,    t=1			4.  r(t)=⟨sin 2t,cos 4t,sin 6t⟩,   t=π/2

5. Compute s(t)=(∫_0^t)||r'(u)||du for r(t)=⟨(t^2),2(t^2),(t^3)⟩ and interpret Mathematica's result.

6. For r(t)=⟨4t,1-3t,(Sqrt[24])t⟩, compute s(t) as in the previous exercise.  Then use s(t) to find an arc length parametrization of r(t), that is, find φ(s)=t, where φ is the inverse of s(t), and check to see that r(φ(s)) has unit speed, that is, ||r'(φ(s))||=1.  Lastly, plot r(t) and r(φ(s)) and compare them.

7.  Consider the helix r(t)=⟨a sin t,a cos t,c t⟩. 
a. Find a formula for the arc length of one revolution of r(t).
b. Suppose a helix has radius 10, height 5, and makes three revolutions.  What is its arc length?

8. The Cornu spiral is defined by r(t)=⟨x(t),y(t)⟩, where x(t)=(∫_0^t)sin(((u^2)/2))du and y(t)=(∫_0^t)cos(((u^2)/2))du.
a.  Plot the Cornu spiral over various intervals for t.
b.  Find a formula for its arc length along the interval -a≤t≤a, where a is a positive real number.
c.  What is its arc length in the limit as a->∞?

13.4.  Curvature

Students should read Section 13.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Vector tools previously studied including arc length enables one to study the idea of curvature, which serves as a measure of how a curve bends, that is, the rate of change in direction of a curve.  In arriving at a definition of curvature, consider a path in vector form and parametrized by

r(t) = ⟨x(t), y(t), z(t)⟩

The parametrization is classified as regular if r' (t) ≠  0 for all  values of t and for which r (t) is defined.  Assume then that r(t) is regular and define the unit tangent vector in the direction of r' (t), denoted T (t), as follows:

T(t)=((r'(t))/(||r'(t)||)).

This unit tangent vector T at any point enables us to determine the direction of the curve at that point, so one may define the curvature κ (Greek letter kappa) at a point as

κ=||((dT)/ds)||=((||T'(t)||)/(||r'(t)||)),

which represents the magnitude of the rate of change in the unit tangent vector with respect to arc length. One denotes the vector dT/ds as the curvature vector.  Its scalar length therefore measures curvature.  For example, a straight line has κ=0 (zero curvature) as one would expect.  For a circle of radius ρ, we have κ=1/ρ (reciprocal of ρ).  This makes sense since a larger circle should have smaller curvature.  In general, if we were to secure a circle, called the osculating circle, that best fits a curve at a specific point on the curve, then curvature of the curve at such a point should agree with the curvature of the osculating circle, that is,

κ=(1/ρ)

Moreover, the radius ρ of this circle is called the radius of curvature. Note that the equations linking κ and ρ illustrate their inverse relationship:

κ=(1/ρ) and ρ=(1/κ)

Example 13.9.  Compute the curvature κ for a circle of radius ρ defined by

r(t)=⟨ρ cos t,ρ sin t⟩

Solution: We first compute the unit tangent vector T using the formula T(t)=((r'(t))/(||r'(t)||)):

Clear[r,T,t,ρ]

r[t_]={ρ Cos[t],ρ Sin[t]}
r'[t]
T[t_]=r'[t]/Simplify[Norm[r'[t]]]

Observe that in this output Mathematica is not able to reduce the expression inside the radical, which simplifies to ρ as a result of the fundamental trigonometric identity (cos^2)x+(sin^2)x=1.  This is due to the Norm command, which employs absolute values.  To remedy this, we use the formula ‖r'(t)‖=(Sqrt[r'(t)·r'(t)]) instead of the Norm command.

T[t_]=r'[t]/Sqrt[Simplify[r'[t].r'[t]]]

We then compute the curvature using the formula κ=((||T'(t)||)/(||r'(t)||)):

κ=Sqrt[Simplify[T'[t].T'[t]]/Simplify[r'[t].r'[t]]]

Since the radius ρ is assumed to be positive, we conclude that κ=(Sqrt[(1/(ρ^2))])=|(1/ρ)|=(1/ρ) as expected.

Example 13.10.  Compute the curvature κ for the curve defined by f(x)=(x^2) at the point (3,9).

Solution: Observe that the graph of a function y=f(x) can be parametrized by x=t and y=f(t) and hence r(t)=⟨t,f(t)⟩. In this case the formula for curvature reduces to

Clear[r,t,f]
r[t_]={t,f[t]}

T[t_]=((r'[t])/(Sqrt[r'[t].r'[t]]))
κ= Sqrt[Simplify[T'[t].T'[t]]/Simplify[r'[t].r'[t]]]

which is the same as κ=((|f''(x)|)/(((1+(((f'(x)))^2)))^(3/2))). With f(x)=(x^2), we get

f[t_]=(t^2)
κ

At x=t=3, the curvature becomes

κ/.t->3

Here is a plot of the curvature along with the function.

Plot[{f[t],κ}, {t,0,3}]

Example 13.10.  Compute the curvature κ and the radius of curvature ρ for the curve defined by

r(t)=⟨1-t,(t^2)+1, (2/3)(t^3)+1⟩ at t=1/2.

Solution: Again we begin by computing the unit tangent vector T:

Clear[r,T,t,κ]

r[t_]={1-t,t^2+1,(2/3)t^3+1}
r'[t]
T[t_]=r'[t]/Sqrt[Simplify[r'[t].r'[t]]]

We then compute the curvature using the same formula as in the previous example and evaluate it at t=1/2:

κ=Sqrt[Simplify[T'[t].T'[t]]/Simplify[r'[t].r'[t]]]
κ/.t->1/2

Hence, the curvature κ=8/9 at t=1/2 and the corresponding radius of curvature is ρ=1/κ=9/8.

Curvature Formula (Cross Product)

There is an alternative formula for calculating the curvature of space curves that involves the cross product and eliminates the need to compute the unit tangent vector function:

κ=((|| r''(t) × r'(t) ||)/(|| r'(t) ((||)^3)))=((|| a(t) × v(t) ||)/(|| v(t) ((||)^3)))

Example 13.11.  Compute the curvature κ(t) and the radius of curvature for the helix defined by r(t)=⟨cos t,sin t,t⟩ for any real number t.

Solution: We first find the derivative of the unit tangent vector with respect to t.

Clear[r,T,t,κ]
r[t_]={Cos[t], Sin[t],t}
r'[t]
r''[t]

κ[t_]=Sqrt[Simplify[Cross[r''[t],r'[t]].Cross[r''[t],r'[t]]]]/((Sqrt[Simplify[r'[t].r'[t]]])^3)

It follows that κ=(1/2) and ρ=2 for all values of t.  Hence, our helix is a curve of constant curvature.

Exercises

In Exercises 1 and 2, find r'(t) and T(t) and evaluate T(2).
1.  r(t)=(3+2t)i+(2-5t)j+9t k			2. v(t)=⟨sin t,cos t,1⟩

3. Use Mathematica to find the curvature function k(x) for y=cos x. Also plot κ(x) for 0≤x≤1. Where does the curvature assume its maximum value?

4. Determine the unit normal vectors to r(t)=[t i+sin t j] at t=(π/4) and t=((3π)/4).

5. Determine the curvature of the vector-valued function r(t)=(3+2t)i+6t j+(5-t)k.

6. Find a formula for the curvature of the general helix r(t)=a cos t i+a sin t j+c t k.

13.5.  Motion in Three Space

Students should read Section 13.5 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Recall that the velocity vector is the rate of the change of the position vector with respect to time while the acceleration vector represents the rate of change of the velocity vector with respect to time.  Moreover, speed is defined to be the absolute value of the velocity vector. In short, we have the following:

v(t)=r'(t), s(t)=‖v(t)‖ and a(t)=v'(t)=r''(t)

One can secure the velocity vector and the position function if the acceleration vector is known via integration.  More specifically:

v(t)=(∫_0^t)a(u)du+(Subscript[v, 0]) where (Subscript[v, 0]) represents the initial velocity vector and r(t)=(∫_0^t)v(u)du+(Subscript[v, 0])t+(Subscript[r, 0]) where (Subscript[r, 0]) is the initial position.

Example 13.12.  Find the velocity vector, the speed, and the acceleration vector for the vector-valued function r(t)=(t^3)i+(1-t)j+4(t^2)k at time t=1.

Solution:

Clear[r,v,s,a]
r[t_]={t^3,1-t,4t^2}
v[t_]=r'[t]
s[t_]=Sqrt[v[t].v[t]]
a[t_]=r''[t]
v[1]
s[1]
a[1]

Thus, v(1)=r'(1)=3i-j+8k, s(1)=(Sqrt[74]), and a(1)=6i+8k.

Example 13.13.  Find r(t) and v(t) if a(t)=t i+4j subject to the initial conditions v(0)=3i-2j and r(0)=0.

Solution: We first solve for v(t) by integrating a(t):

Clear[r,v,a]
a[t_]={t,4}
v[t_]=Integrate[a[u],{u,0,t}]+{v01,v02}

Here, the constant vector of integration (Subscript[v, 0])=⟨(Subscript[v, 01]),(Subscript[v, 02])⟩=⟨3,-2⟩ equals the initial velocity:

Solve[v[0]=={3,-2},{v01,v02}]

Thus, v(t)=((t^2)/2)i+4t j+(3i-2j).

v[t_]=v[t]/.{v01->3,v02->-2}

Next, we solve for r(t) by integrating v(t):

r[t_]=Integrate[v[u],{u,0,t}]+{r01,r02}

Again, the constant vector of integration (Subscript[r, 0])=⟨(Subscript[r, 01]),(Subscript[r, 02])⟩=⟨0,0⟩ equals the initial position:

Solve[r[0]=={0,0},{r01,r02}]

Hence, r(t)=(((t^3)/6)+3t)i+(2(t^2)-2t)j.

Components of Acceleration

There are two components of acceleration: tangential and normal.  More precisely, the acceleration vector a can be decomposed as a=(Subscript[a, T])T+(Subscript[a, N])N, where (Subscript[a, T])=(((d^2)s)/(dt^2))=((a · v)/(|| v ||)) is the tangential component and (Subscript[a, N])=((κ((ds/dt)))^2)=((|| a × v ||)/(|| v ||)) is the normal component.  Moreover, one has (Subscript[a, T]^2)+(Subscript[a, N]^2)=((‖a‖)^2) so that (Subscript[a, N])=(Sqrt[||a((||)^2)-((Subscript[a, T])^2)]) and (Subscript[a, T])=(Sqrt[||a((||)^2)-((Subscript[a, N])^2)]).

Example 13.14.  Determine the tangential and normal components of acceleration for the vector function r(t)=⟨(t^3),(t^2),t⟩.

Solution:

Clear[r,v,s]
r[t_]={t^3,t^2,t}
r'[t]
r''[t]

speed=Simplify[Sqrt[r'[t].r'[t]]]

The result in the last output represents the speed at time t.  In order to secure the tangential component of the acceleration, we differentiate the previous output:

at=D[speed,t]

The normal component of the acceleration is

an=(Sqrt[r''[t].r''[t]-(at^2)])

Simplify[an]

NOTE: The components of acceleration can also be found through the formulas (Subscript[a, T])=((a · v)/(|| v ||)) and (Subscript[a, N])=((|| a × v ||)/(|| v ||)), confirmed using Mathematica as follows:

at=r''[t].r'[t]/Sqrt[r'[t].r'[t]]
an=Sqrt[Cross[r''[t],r'[t]].Cross[r''[t],r'[t]]]/Sqrt[r'[t].r'[t]]

Exercises

In Exercises 1 and 2, calculate the velocity and acceleration vectors and the speed at the time indicated:
1.  r(t)=(t^2)i+(1-t)j+(5(t^2))k, t=2.		2.  r(t)=cos t i+sin t j+tan (2t)k,t=(π/6).

3.  Sketch the path r(t)=(1-(t^2))i+(1-t)j for -3≤t≤3 and compute the velocity and acceleration vectors at t=0, t=1, and t=2.

4.  Find v(t) given a(t) and the initial velocity (Subscript[v, 0]).
a.  a(t)=t i+3j,v(0)=(1/2)i+2j 			b.  a(t)=(e^(2t))i+0j+(t+2)k,v(0)=i-3j+2k

5.  Find r(t) and v(t) given a(t) together with the initial velocity and position at rest:
a.  a(t)=(e^(3t))i+4t j+(t-2)k,v(0)=0i+0j+k,r(0)=0i+3j+4k
b. a(t)=0i+0j+sin t k,v(0)=i+j,r(0)=i.

6.  Find the decomposition of a(t) into its tangential and normal components at the indicated point:
a.  r(t)=(3-4t)i+(t+1)j+(t^2)k at t=1
b.  r(t)=t i+(e^(-t))j+t (e^(-t))k at t=0

7.  Show that the tangential and normal components of acceleration of the helix given by r(t)=(cos t) i+(sin t) j +t k are equal to 0 and 1, respectively.

Chapter 14	 Differentiation in Several Variables

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

14.1  Functions of Two or More Variables

Students should read Section 14.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

14.1.1 Plotting Level Curves using ContourPlot

We begin with plotting level curves f(x,y)=c of a function of two variables. The command to plot level curves is ContourPlot[f,{x,a,b},{y,c,d}].

Most of the options for ContourPlot are the same as those for Plot.  In the following example, we consider the option ImageSize.

Example 14.1.  Plot the level curves of f(x,y)=(x^2)+x y -(y^2) .

Solution:  Let us first plot the level curves using the default settings of Mathematica.

Clear[x,y,f]
f[x_,y_]:=(x^2)+x y -(y^2)

ContourPlot[f[x,y],{x,-5,5},{y,-5,5},ImageSize->{250}]

To get the level curves on the xy-plane without the shading, the colors, and the frame, but with the coordinate axes, we use the following options of ContourPlot.

ContourPlot[f[x,y],{x,-5,5},{y,-5,5},Frame->False,Axes->True,ContourShading->False,ImageSize->{250}]

Contours is an option of ContourPlot that can be used in two different ways: Contour->n displays n equally spaced contour curves while Contour->list plots level curves f(x,y)=c where c is an element of the list list.

To plot 15 level curves, we evaluate

ContourPlot[f[x,y],{x,-1,1},{y,-1,1}, Contours->15,ImageSize->{250} ]

Here is an example when list={-10,-5,-2, -1, 0,1, 2,  5, 10}.

ContourPlot[f[x,y],{x,-5,5},{y,-5,5},Contours->{-10,-5,-2,-1, 0, 1, 2, 5, 10},ImageSize->{250}]

14.1.2  Plotting Surfaces using Plot3D

Plot3D is the three-dimensional analog of the Plot command. Given an expression in two variables and the domain for the variables, Plot3D produces a surface plot.

The basic syntax to plot the graph of a function of two variables is Plot3D[ f,{x, a, b},{y, c, d}], where f is a function of x and y  with  a≤x≤b and c≤y≤d.

The command to plot the graphs of two or more functions on the same coordinate axes is Plot3D[{f, g, h, .... }, {x, a, b}, {y, c, d}], where f, g, h, ...  are the functions to be plotted.

We will begin with the default settings of plotting a graph of a function of two variables.

Example 14.2.  Plot f(x,y)=sin(x-cos y).

Solution:

Plot3D[Sin[x-  Cos[y]],{x,-3,3},{y,-3,3}]

Example 14.3.  Plot the graphs of f(x,y)=3x+4y-3 and g(x,y) =10 sin(xy) on the same axes.

Solution:  We will use red color for the graph of f and blue for that of g. This is given using the option PlotStyle.

Plot3D[{3x+4y-3,10Sin[x y]},{x,-5,5},{y,-5,5},PlotStyle->{Green,Blue},ImageSize->{250}]

NOTE: One of the most significant improvements of Mathematica 7.0 over the previous editions is its graphics capability.  Plot3D has many options.  Evaluate the command Options[Plot3D] to see the many options you have to plot a nice graph.

We will discuss some of these options below.

ViewPoint

In Mathematica 7.0, we can rotate the graph of a function of two variables by simply clicking on the graph and dragging the mouse around to see the graph from any point of view.

The option ViewPoint specifies the point in space from which the observer looks at a graphics object. The syntax for choosing a view point of a surface is Plot3D[f[x, y], {x, a, b}, {y, c, d}, ViewPoint->{A, B, C} ]. The default value for {A, B, C}  is {1.3,-2.4,2.0}. This may be changed by entering values directly.

To view a graph from directly in front {0, -2,0}; in front and up {0, -2,2}; in front and down {0, -2,-2}; left hand corner {  -2,-2,0};  directly above  {0, 0,2}.

Plot3D[ f[x, y], {x, a, b}, {y, c, d}, ViewPoint -> view ] produces a plot viewed from view. The possible values of view are Above (along positive z-axis), Below (along negative z-axis), Front (along negative y-axis), Back (along  positive y-axis), Left  (along the negative x-axis), and Right (along the positive x-axis).

Example 14.4.  Plot f(x,y)=cos x sin y using ViewPoint option to view the graph from various view points.

Solution: We leave it to the reader to navigate all of the above choices. We will consider a few of them.

Clear[f]
f[x_,y_]=Cos[x]  Sin[y]

Here is a plot of the graph using the default setting for ViewPoint:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi} ,PlotRange->All,
ImageSize->{250} ]

View from directly in front:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi},ViewPoint->Front ,
PlotRange->All ,ImageSize->{250}]

View from in front and up:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi},ViewPoint->{0,-2,2},
PlotRange->All,ImageSize->{250}]

View from in front and down:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi},ViewPoint->{0,-2,-2},
PlotRange->All,ImageSize->{250}]

View from directly above:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi},ViewPoint->Above,
PlotRange->All,Ticks->{Automatic,Automatic,{-1,0,1}},
ImageSize->{250}]

View from the right:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi},ViewPoint->Right,
PlotRange->All,ImageSize->{250}]

NOTE: As we pointed out earlier, we can also select different viewpoints by clicking on the graph and dragging the mouse around until we get the desired viewpoint.

Mesh, MeshStyle, MeshShading

The option Mesh specifies the type of mesh that should be drawn.

The option MeshStyle specifies the style in which a mesh should be drawn.

The option MeshShading is an option for specifying a list of colors to be used between mesh divisions.

We illustrate some uses of these options in the example below.

Example 14.5.  Plot f(x,y)=cos x sin y using various options involving Mesh.

Solution:

Clear[f]
f[x_,y_]=Cos[x]  Sin[y]

To plot a graph without a mesh we use the setting Mesh->None.

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi} , Mesh->None,
ImageSize->{250} ]

Mesh->n plots a surface with only n×n meshes.

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi} , Mesh->8 ,
ImageSize->{250}]

We can choose the color of the mesh using MeshStyle.

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi} , MeshStyle->
{Red,Black},ImageSize->{250} ]

Here is another use of MeshStyle:

Plot3D[ f[x,y],{x,-2Pi,2Pi},{y,-2Pi,2Pi} , MeshStyle->
{Dashing[0.01],None},ImageSize->{250} ]

To display a plot with selected colors between meshes we use MeshShading:

Plot3D[f[x,y],{x,-2Pi, 2Pi},{y,-2Pi,2Pi},MeshShading->{{Blue, Red, White},{Purple,Green,Black}},ImageSize->{250} ]

Here is a neat example in Mathematica 7.0:

Plot3D[(x^2-y^2)/(x^2+y^2)^2,{x,-1.5,1.5},{y,-1.5,1.5},BoxRatios->Automatic,PlotPoints->25,MeshFunctions->{#3&},MeshStyle->Purple,MeshShading->{None,Green,None,Yellow},ImageSize->{250}]

BoxRatios

The option BoxRatios specifies the ratio of the lengths of the sides of the box. This is analogous to specifying the AspectRatio of a two-dimensional plot.  For Plot3D, the default setting is BoxRatios->Automatic.

Example 14.6.  Plot f(x,y)=(e^(1-(x^2)-(y^2))) using the BoxRatio option.

Solution:

Clear[f]
f[x_,y_]= (E^(1-(x^2)- (y^2)))

Plot3D[ f[x,y],{x,-2,2},{y,-2,2},ImageSize->{250}]

Plot3D[ f[x,y],{x,-2,2},{y,-2,2},BoxRatios->{1,1,0.62`},
ImageSize->{250}]

AxesLabel

The option AxesLabel  is a command used to label the axes in plotting.

Example 14.7.  Plot f(x,y)=(Sqrt[9-(x^2)-(y^2)]) using the AxesLabel option.

Solution:

Clear[f]
f[x_,y_]= (Sqrt[9-(x^2)-(y^2)])

Plot3D[f[x,y],{x,-3,3},{y,-3,3},AxesLabel->{"x ","y ","z "},
ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

NOTE: To label a graph, use the PlotLabel option as shown following:

Plot3D[f[x,y],{x,-3,3},{y,-3,3}, AxesLabel->{"x ","y ","z "},
PlotLabel->"Upper hemisphere",BoxRatios->Automatic,ImageSize->{250},
ImagePadding->{{15,15},{15,25}}]

ColorFunction

The option ColorFunction specifies a function to apply to the values of the function being plotted to determine the color to use for a particular region on the xy-plane. It is an option for Plot3D, ListPlot3D, DensityPlot, and ContourPlot.  The default setting for ColorFunction is ColorFunction->Automatic.  ColorFunction->Hue yields a range of colors.

Example 14.8.  Plot f(x,y)=sin((x^2)+(y^2)) +(e^(1- ((x )^2)-(y^2)))  in various colors using the ColorFunction option.

Solution:

Clear[f,x,y]
f[x_,y_]= Sin[(x^2)+(y^2)]+ (E^(1-(x^2)- (y^2)))

Plot3D[f[x,y],{x,-Pi,Pi},{y,-Pi,Pi},ColorFunction->Hue,
ImageSize->{250}]

Here are other ways to use ColorFunction.

Plot3D[f[x,y],{x,-Pi,Pi},{y,-Pi,Pi},ColorFunction-> "Rainbow",
ImageSize->{250} ]

Plot3D[f[x,y],{x,-Pi,Pi},{y,-Pi,Pi},ColorFunction->
"BlueGreenYellow",ImageSize->{250}]

NOTE: We can use PlotStyle option to select color for graphs. The plot below uses this option.

Plot3D[f[x,y],{x,-Pi,Pi},{y,-Pi,Pi},PlotStyle-> Yellow,
ImageSize->{250} ]

RegionFunction

The option RegionFunction specifies the region to include in the plot drawn.

Example 14.9.  Plot f(x,y)={({{10sin(3x-y),, if  (x^2)+(y^2)<4;}, {(x^2)+(y^2)-5,    , otherwise}}).

Solution: We will use the command RegionFunction to specify the domain (x^2)+(y^2)<4 as follows.  Note that we have used Show to display the graphs.

Clear[plot1,plot2] 
plot1 =Plot3D[10Sin[3x-y],  {x,-4,4},{y,-4,4},PlotStyle->Blue,
 RegionFunction->Function[{x,y,z}, x^2+y^2<4]];
plot2 =Plot3D[  (x^2)+(y^2) -5,{x,-4,4},{y,-4,4},PlotStyle->Red, 
 RegionFunction->Function[{x,y,z}, x^2+y^2 ≥4]];
Show[plot1,plot2,ImageSize->{250}]

If we want to focus on a particular part of a surface defined by a function, we can use the option RegionFunction. The following example shows this point.

Example 14.10.  Plot the graph of f(x,y) =(x^2)-3x y-2(y^2)  and show the portion of the surface direclty above the unit circle centered at the origin.

Solution: We will use the option ViewPoint.

Clear[plot1,plot2,f, x,y]
f[x_,y_]=(x^2) -3 x y -2(y^2) 
plot1 =Plot3D[f[x,y],  {x,-4,4},{y,-4,4},PlotStyle->Blue,
 RegionFunction->Function[{x,y,z},x^2+y^2<1 ]];
plot2 =Plot3D[f[x,y] ,{x,-4,4},{y,-4,4},PlotStyle->Red,
RegionFunction->Function[{x,y,z}, x^2+y^2>1 ]]; 
Show[plot1,plot2 ,ViewPoint->Front,ImageSize->{250} ]

14.1.3  Plotting Parametric Surfaces using ParametricPlot3D

ParametricPlot3D is a direct analog of ParametricPlot.  Depending on the input, ParametricPlot3D  produces a space curve or a surface. ParametricPlot3D[{f, g, h}, {t, a, b }] produces a three-dimensional space curve parametrized by the variable t, which runs from  a to b. ParametricPlot3D[{f, g, h}, {t, a, b },{u, c, d}] produces a two-dimensional surface parametrized by t and u. Options are given to ParametricPlot3D the same way as for  Plot3D. Most of the options are the same.

Example 14.11.  Plot the curve that is parametrized by x=sin t, y=cos t and z=t/3 with 0≤t≤2π.

Solution:

ParametricPlot3D[{Sin[t],Cos[t],(t/3)},{t,0,2 π},ImageSize->{250},
ImagePadding->{{15,15},{15,15}}]

Example 14.12.  Plot the surface that is parametrized by  x = u cos u (4 + cos(u + v)), y = u sin u (4 + cos(u + v)), and z = u sin(u + v).

Solution:

ParametricPlot3D[{u Cos[u] (4+Cos[u+v]),u Sin[u] (4+Cos[u+v]),u Sin[u+v]},{u,0,4 π},{v,0,2 π},ImageSize->{250}]

14.1.4  Plotting Level Surfaces using ContourPlot3D

ContourPlot3D is the command used to plot level surfaces of functions of three variables.  Its syntax is ContourPlot3D[f,{x,a,b}, {y,c,d},{z,e,f}].  Most of the Options for ContourPlot3D are the same as those of Plot3D.  Below we will consider the option Contours of ContourPlot3D.

Example 14.13.  Plot level surfaces of f(x,y,z)=(x^2)+(y^2)+(z^2).

Clear[x,y,z,f]
f[x_,y_,z_]=(x^2)+(y^2)+(z^2)
ContourPlot3D[f[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},ImageSize->{250}]

The following displays five (5) equally spaced contour surfaces of f.

ContourPlot3D[f[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},
Contours->5,ImageSize->{250}]

The following  displays three level surfaces f(x,y,z)=c, where c=1, 4, 9.

ContourPlot3D[f[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},
Contours->{1,4,9},ImageSize->{250}]

Notice that we only see one sphere. The other two are enclosed in the sphere of radius 3 corresponding to c=9.  One way to remedy this is to plot the level surfaces one by one. For this we use the GraphicsArray command. First, let us define the level surfaces as function of c:

Clear[c,plot]
plot[c_]:=ContourPlot3D[f[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},
Contours->{c}]

Here are the three level surfaces corresponding to c=1, 4, 9.

Show[GraphicsArray[{plot[1],plot[4],plot[9]}]]

Exercises

In Exercises 1 through 4, plot the level curves and the graphs of the given functions.
1.  f(x,y)=x (y^5)-(x^5)y for -10≤x≤10, -10≤y≤10
2.  f(x,y)=(((x^2)+2y)/(1+(x^2)+(y^2))) for -10≤x≤10, -10≤y≤10
3.  f(x,y)=(sin y )(e^(cos x)) for -2π≤x≤2π, -2π≤y≤2π
4.  f(x,y)= sin(x+sin(y)) for -4π≤x≤4π, -4π≤y≤4π

In Exercises 5 through 7, use at least two nondefault options to plot the given functions.
5.  f(x,y)=sin (x-2y)(e^(1/(y-x))) for  -2π≤x≤2π, -2π≤y≤2π
6.  f(x,y)=4-3|x|-2|y| for -10≤x≤10, -10≤y≤10
7.  f(x,y)=(tanh^(-1))(x/y) for -5≤x≤5, -5≤y≤5

8.  Plot f(x,y)={({{(x^2)+(y^2)-4      if (x^2)+(y^2)<4}, {4-(x^2)+3(y^2)         otherwise }})

9.  Plot the portion of the helicoid (spiral ramp) that is defined by:
x = u cos v,  y = u sin v, z = v for 0 ≤ u ≤ 3  and  -2π ≤ v ≤2 π

10. Use ContourPlot3D to plot the level surfaces of the function f(x,y,z)=9-(x^2)-(y^2)-(z^2).

14.2  Limits and Continuity

Students should read Section 14.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

14.2.1  Limits

If f(x,y) is a function of x and y, and if the domain of f contains a circle around the point (a,b), we say that the limit of f at (a,b) is L if and only if f(x,y) can be arbitrarily close to L for all (x,y) arbitrarily close (a,b).

More precisely, for a given ϵ>0, there exists a δ>0 such that for every (x,y) is in the domain of f,

0<(Sqrt[(((x-a))^2)+(((y-b))^2)])<δ  ⟹ |f(x,y)-L|<ϵ

If this is the case, we write

(lim_((x,y)->(a,b)))f(x,y)=L

The Limit command of Mathematica is restricted to functions of one variable. However, we can use it twice to find the limit of function of two variables provided the limit exists.

Example 14.14.  Find (lim_((x,y)->(3,4)))((x^2)+(y^2)).

Solution: We can easliy determine that the limit exists. We can find the limit by evaluating

Limit[Limit[(x^2)+(y^2),x->3],y->4]

The plot following confirms this.

Clear[plot1,plot2]
plot1=Plot3D[(x^2)+(y^2),{x, 1,4},{y, 3, 5}];
plot2=Graphics3D[{Red,PointSize[.025],Point[{3,4,25}]}];
Show[plot1,plot2,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

Example 14.15.   Find (lim_((x,y)->(4,1)))((3x+(y^2))/(x-4y)).

Solution: We will evaluate the limit in two different orders. The limit in which we use limit with x first and then with y is

Clear[f,x,y]
f[x_,y_]=((3x+(y^2))/(x-4y))

The limit in which we use limit with x first and then with y is

Limit[Limit[f[x,y],x->4],y->1]

The limit in which we use limit with y first and then with x is

Limit[Limit[f[x,y],y->1],x->4]

Here is the plot of the graph near the point (4,1). Observe that the graph of the function is in green and the point (4,1,0) is in red. For a better comaprison, we have colored the xy-plane light blue. You may need to rotate the graph to see the point (4,1,0) on the xy-plane and see how the graph behaves when (x,y) is close to (4,1).

Clear[plot1,plot2]
plot1=Plot3D[{f[x,y],0},{x,3,5},{y, 0, 2},PlotStyle->{Green, LightBlue},PlotPoints->100];
plot2=Graphics3D[{Red,PointSize[.025],Point[{4,1,0}]}];
Show[plot1,plot2,ImageSize->{250},
ImagePadding->{{15,15},{15,15}}]

Here is the animation with x as the animation parameter.

Important Note: If you are reading the printed version of this publication, then you will not be able to view any of the animations generated from the Animate command in this chapter.  If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then evaluate each Animate command to view the corresponding animation.  Just click on the arrow button to start the animation.  To control the animation just click at various points on the sliding bar or else manually drag the bar.

Animate[Plot[f[x,y],{y, 0,3},PlotRange->{-20,20}],{x,3, 5}]

Example 14.16.  Find (lim_((x,y)->(0,0)))((sin x sin y)/(x y)).

Solution: We will evaluate the limit in two different orders.

Clear[f,x,y]
f[x_,y_]=((Sin[x  y])/(x y))

Limit[Limit[f[x,y],x->0],y->0]

Limit[Limit[f[x,y],y->0],x->0]

Here is the plot of the graph and the point (0,0,1).

Clear[plot1,plot2]
plot1=Plot3D[f[x,y],{x,-1,1},{y, -1, 1},PlotStyle-> Green];
plot2=Graphics3D[{Red,PointSize[.02],Point[{0,0,1}]}]; 
Show[plot1,plot2,ImageSize->{250}]

If we rotate this graph to a suitable position, we notice that the limit exists.  Here are animations with x and y as animation parameters, respectively.

Animate[Plot[f[x,y], {x, -2,2},PlotRange->{0,1}], {y, -2, 2}]

Animate[Plot[f[x,y], {y, -2,2},PlotRange->{0,1}], {x, -2, 2}]

Example 14.17.  Find (lim_((x,y)->(0,0))) x ln y.

Solution:

Clear[f,x,y]
f[x_,y_]=x Log[y]

Limit[Limit[f[x,y],x->0], y->0]

Limit[Limit[f[x,y],y->0], x->0]

Clear[plot1,plot2]
plot1=Plot3D[{f[x,y],0},{x,-1,1},{y, -1, 1},PlotStyle->{ Green,LightBlue}];
plot2=Graphics3D[{Red,PointSize[.025],Point[{0,0,0}]}];
Show[plot1,plot2,ImageSize->{250},
ImagePadding->{{15,15},{15,15}}]

Here is the animation with x as the animation parameter.

Animate[Plot[f[x,y], {y, -2,2},PlotRange->{-10,10}], {x, -2, 2}]

Example 14.18.  Consider the function f(x,y)=((x (y^2))/((x^2)+(y^4))). Show that (lim_((x,y)->(0,0)))f(x,y) does not exist.

Solution:

Clear[f,x,y]
f[x_,y_]=((x (y^2))/((x^2)+(y^4)))

Limit[Limit[f[x,y],x->0], y->0]

Limit[Limit[f[x,y],y->0], x->0]

Limit[Limit[f[x,y],y-> m x], x->0]

However, note that the limit along the curve y=(Sqrt[x]) is

Limit[Limit[f[x,y],y->  (Sqrt[x])], x->0]

Hence, the limit does not exist.  Here is the plot of the function:

Plot3D[f[x,y],{x,-1,1},{y, -1, 1},ImageSize->{250}]

14.2.2  Continiuty

Recall that a function f of two variables x and y is continuous at the point (a,b) if and only if (Subscript[limi, (x,y)->(a,b)])f(x,y)=f(a,b).

Example 14.19.  Let f(x,y)={({{1-(x^2)-(y^2), if (x^2)+(y^2)<1}, {0,                 if (x^2)+(y^2)≥1}}).  Is f continuous?

Solution: Clearly, f is continuous at all points inside and outside the circle of radius 1.  To check continuity on the unit circle, we let x=r cos t and y=r sin t. We then let r->1.

Clear[x,y,r,s,t,f] 
f[x_,y_]=1-(Sqrt[(x^2)+(y^2)])

x=r Cos[t]
y=r Sin[t]

Simplify[f[x,y]]

Limit[f[x,y],r->1]

The command below evaluates f on the circle.

Simplify[f[x,y]/.r->1]

Thus, the limit and the value of f are equal at all points on the unit circle. Hence, f is continuous everywhere.  Here is the graph.

Clear[plot1,plot2]
plot1=Plot3D[f[x,y],{x,-5,5},{y,-5,5},PlotStyle->Red, 
 RegionFunction->Function[{x,y,z},x^2+y^2<1],Mesh->None];
plot2=Plot3D[  0,{x,-5,5},{y,-5,5},PlotStyle->LightBlue,
 RegionFunction->Function[{x,y,z}, x^2+y^2 ≥1],Mesh->None];
Show[plot1,plot2,ImageSize->{250}]

Exercises

In Exercises 1 through 4, find the limit, if it exists.
1.  (lim_((x,y)->(1,-1)))(2(x^2)y +x (y^2))		2. (lim_((x,y)->(1,1)))((3(x^2)+(y^2))/((x^2)-y))
3.   (lim_((x,y)->(0,0)))((tan x sin y)/(x y))	 		4.  (lim_((x,y)->(0,0))) sin x ln y

5. Consider the function f(x,y)=(((x^2) +(y^2))/((x^2)+(y^4))). Show that (lim_((x,y)->(0,0)))f(x,y) does not exist.

6.  Let f(x,y)={({{(x^2)-(y^2), if x+y<0}, {2x+y,  if x+y≥0}}).
Is f continuous?

7.  Let f(x,y)=((x y)/((x^2)+(y^2))). The domain of f is the whole plane without the origin. Is it possible to define f(0,0) so that f is continuous everywhere? Plot the graph of f to support your conclusions.

8.  The domain of f(x,y)=((x y)/(x+y))  is the whole plane without the line y=-x.   Is it possible to define f(0,0) so that f is continuous everywhere? Plot the graph of f to support your conclusions.

14.3  Partial Derivatives

Students should read Section 14.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Recall that the Mathematica command for the partial derivative of a function f with respect to x is D[f, x], and D[f,{x,n}] gives the nth partial derivative of f with respect to x. The multiple (mixed) partial derivative of f with respect to (Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),... is obtained by D[f,(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),...].  We can access this command from  BasicMathInput. The symbols are (∂_□) □ and (∂_(□,□))□ .

Example 14.20.  Find the first partial derivatives of (x^3)+(y^2) with respect to x and (y^( )).

Solution:  We give two methods of input.

Method 1: We can type all the inputs and the command as follows:

Clear[x,y]
D[x^3+y^2,x]

D[x^3+y^2,y]

Method 2:  We can use the BasicInput palette to enter the inputs.

(∂_x) ((x^3)+(y^2))

(∂_y) ((x^3)+(y^2))

Example 14.21.  Find the four second partial derivatives of (x^3)sin(y)+(e^(x y)).

Solution: Let z=(x^3)sin y+(e^(x y)).  We again demonstrate two methods of input.

Method 1:

We can find (Subscript[z, xx]) by

Clear[x,y]
D[x^3*Sin[y]+E^(x*y),{x,2}]

We can find (Subscript[z, yy]) by

D[x^3*Sin[y]+E^(x*y),{y,2}]

We can find (Subscript[z, xy]) by

D[x^3*Sin[y]+E^(x*y),x,y]

(Subscript[z, yx])  is given by

D[x^3*Sin[y]+E^(x*y),y,x]

NOTE: Clairaut's Theorem states that if the mixed partial derivatives (Subscript[f, x y]) and (Subscript[f, y x]) are continuous at a point (x,y), then they are equal: (Subscript[f, x y])=(Subscript[f, y x]).  The last two outputs confirm Clairaut's Theorem for this particular example.

Method 2: Here is the input using the palette symbol (∂_(□,□))□:

Clear[x,y]
(∂_(x,x))((x^3)*Sin[y]+(E^(x*y)))   
(∂_(y,y))((x^3)*Sin[y]+(E^(x*y)))
(∂_(x,y))((x^3)*Sin[y]+(E^(x*y)))
(∂_(y,x))((x^3)*Sin[y]+(E^(x*y)))

Example 14.22.  Evaluate the first partial derivatives of x y+y (z^2)+x z at (-1,2,3).

Solution: Recall that Expr /.{(Subscript[x, 1])->(Subscript[a, 1,])(Subscript[x, 2])->(Subscript[a, 2,])(Subscript[x, 3])->(Subscript[a, 3,])... } is the command for substituting (Subscript[x, 1]) by (Subscript[a, 1]), (Subscript[x, 2]) by (Subscript[a, 2]), (Subscript[x, 3]) by (Subscript[a, 3]), .... , in Expr.

Clear[x,y,z]
D[x*y + y*z^2 + x*z,x]/.{x-> -1, y->2, z->3}

D[x*y + y*z^2 + x*z,y]/.{x->-1, y->2, z->3}

D[x*y+y*z^2+x*z,z]/.{x->-1, y->2, z->3}

Example 14.23.  Let f(x,y,z)=y (e^x)+x (e^(-y))ln z. Find   (Subscript[f, x x x]), (Subscript[f, x y z]), (Subscript[f, x z z]), (Subscript[f, z x z ]), and (Subscript[f, z z x]).

Solution: First, we define f(x,y,z) in Mathematica. We can use the (∂_(□,□))□ notation. Since the palette gives only two boxes for the variables, we need to add one more box. This can be done by using CTRL +, (comma), that is, hold the CONTROL key and press the COMMA button.  Note also that the command D[f[x,y,z],x,y,z] gives (Subscript[f, xyz]). We demonstrate both methods.

Clear[x,y,z,f]
f[x_,y_,z_]:=y*(E^x)+x *Log[z]*(E^(-y))

(∂_(x,x,x)) f[x,y,z]

(∂_(x,y,z)) f[x,y,z]

(∂_(x,z,z)) f[x,y,z]

D[f[x,y,z],z,x,z]

D[f[x,y,z],z,z,x]

Example 14.24.  Let f(x, y) = x y (((x^2)-(y^2))/((x^2)+(y^2))) if (x,y)≠(0,0) and f(0,0)=0. 
a) Find (Subscript[f, x])(x,y) and (Subscript[f, y])(x,y) for (x,y)≠(0,0).
b) Use the limit definition to find (Subscript[f, x])(0,0) and (Subscript[f, y])(0,0).
c) Find (Subscript[f, x y]) (x,y) and (Subscript[f, y x])(x,y) for (x,y)≠(0,0).
d) Use the limit definition to find (Subscript[f, x y])(0,0) and (Subscript[f, y x])(0,0).

Solution: We will first define f using the If command.

Clear[x,y,f,fx,fy,fxy,fyx]
f[x_,y_]=If[{x,y}≠{0,0}, x y (((x^2)-(y^2))/((x^2)+(y^2))),0]

a)  Let fx and fy denote the partial derivatives with respect to x and y, respectively. Then

fx[x_,y_]=D[f[x,y],x]
fy[x_,y_]=D[f[x,y],y]

If we use the FullSimplify command to simplify the preceding output, we get

FullSimplify[fx[x,y]]
FullSimplify[fy[x,y]]

Thus, (Subscript[f, x])(x,y)=((y((x^4)+4(x^3)(y^2)-(y^4)))/((((x^2)+(y^2)))^2)) and  (Subscript[f, y])(x,y)=((x((x^4)-4(x^2)(y^2)-(y^4)))/((((x^2)+(y^2)))^2)) if (x,y)≠(0,0).

b) We use the limit definition (Subscript[f, x])(0,0)=(lim_(h->0 ))((f(0+h,0)-f(0,0))/h)and (Subscript[f, y])(0,0)=(lim_(k->0 ))((f(0,0+k)-f(0,0))/k)  to find the partial derivatives at (0,0).

Clear[h,k]
Limit[((f[0+h,0]-f[0,0])/h),h->0]

Limit[((f[0,0+k]-f[0,0])/k),k->0]

Hence, (Subscript[f, x])(0,0)=0 and (Subscript[f, y])(0,0)=0.

c) To find the mixed second partial derivatives, we use fx and fy from the outputs in part a).  Note that the FullSimplify command is used to to get a simplified form of the mixed partial derivatives.

fxy[x_,y_]=FullSimplify[D[fx[x,y],y]]
fyx[x_,y_]=FullSimplify[D[fy[x,y],x]]

Thus, (Subscript[f, x y])=(((x-y)(x+y)((x^4)+10(x^2)(y^2)+(y^4)))/((((x^2)+(y^2)))^3)) and (Subscript[f, y x])=(((x-y)(x+y)((x^4)+10(x^2)(y^2)+(y^4)))/((((x^2)+(y^2)))^3)) for (x,y)≠(0,0). Note that these two functions are equal for (x,y)≠(0,0) in conformity with Clairaut's Theorem, since both are continuous when (x,y)≠(0,0).

d) We use the limit definition of a partial derivative to compute (Subscript[f, xy])(0,0) and (Subscript[f, yx])(0,0).  Recall that we have defined (Subscript[f, x]) as fx[x,y] and (Subscript[f, y]) as fy[x,y].

Then (Subscript[f, xy])(0,0) is given by

Limit[((fx[0 ,0+k]-fx[0,0])/k),k->0]

and (Subscript[f, yx])(0,0) is given by

Limit[((fy[0+h,0]-fy[0,0])/h),h->0]

Thus, (Subscript[f, xy])(0,0)=-1and (Subscript[f, yx])(0,0)=1.  Note that this implies that the mixed partial derivatives are not continuous at (x,y)=(0,0).  To see this graphically, first consider the following graph of f, which confirms that f has partial derivatives everywhere.

Plot3D[f[x,y],{x,-3,3},{y,-3,3},ImageSize->{250}]

Here are the graphs of (Subscript[f, x]) and (Subscript[f, y]), which now show why the second mixed partials at the origin are not equal.

Clear[plot1,plot2]
plot1= Plot3D[fx[x,y],{x,-3,3},{y,-3,3}, PlotStyle->Red,AxesLabel->{"Graph of z=(Subscript[f, x])", None, None}] ;
plot2= Plot3D[fy[x,y],{x,-3,3},{y,-3,3}, PlotStyle->Blue,AxesLabel->{"Graph of z=(Subscript[f, y])", None, None}] ;
Show[GraphicsArray[{plot1,plot2}],ImageSize->{420}]

In addition, the graphs of (Subscript[f, xy]) and (Subscript[f, yx]) show the mixed partials are not continuous at the origin. This is the main reason why the inequalities of the mixed partials at the origin does not contradict Clairaut's Theorem.

Clear[plot1,plot2]
plot1=Plot3D[fxy[x,y],{x,-3,3},{y,-3,3}, PlotStyle->Red, AxesLabel->{"Graph of z=(Subscript[f, xy])", None, None}];
plot2=Plot3D[fyx[x,y],{x,-3,3},{y,-3,3}, PlotStyle->Blue, AxesLabel->{"Graph of z=(Subscript[f, yx])", None, None}];
Show[GraphicsArray[{plot1,plot2}],ImageSize->{420}]

Exercises

1. Let f(x,y)=((((x-y))^2)/((x^2)+(y^2))).  Find:  
a. (Subscript[f, x]) (1,0)	b. (Subscript[f, y])(1,0)	c. (Subscript[f, xy])		d. (Subscript[f, yx])		e.  (Subscript[f, xxy])

2. Find the first partial derivatives of z=(x^3)(y^2) with respect to (x^( ))and y.

3. Find the four second partial derivatives of (x^2)cos(y) +tan(x (e^y)).

4. Evaluate the first partial derivatives of f(x,y,z)=(e^(-z))xy +(yz^2)+xz at (-1,2,3).

5. Let f(x,y,z)=(((x^4)(y^3))/((z^2)+sin x)).  Find (Subscript[f, xxx]), (Subscript[f, xyz]), (Subscript[f, xzz]), (Subscript[f, zxz ]), and (Subscript[f, zzx]).

6. Let f(x, y) = ((x (y^2))/((x^2)+(y^4))) if (x,y)≠(0,0) and f(0,0)=0. 
a.  Find (Subscript[f, x])(x,y) and (Subscript[f, y])(x,y) for (x,y)≠(0,0).
b.  Use the limit definition to find (Subscript[f, x])(0,0) and (Subscript[f, y])(0,0).
c.  Find (Subscript[f, xy]) (x,y) and (Subscript[f, yx])(x,y) for (x,y)≠(0,0).
d.  Use the limit definition to find (Subscript[f, xy])(0,0) and (Subscript[f, yx])(0,0).

14.4  Tangent Planes

Students should read Section 14.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let z=f(x,y) be a function of two variables.  The equation of the tangent plane at the point (a,b,f(a,b)) is given by

z=(Subscript[f, x])(a,b)(x-a) +(Subscript[f, y])(a,b)(y-b)+f(a,b)

Example 14.25.  Let f(x, y) =  (x^2)+(y^2).
a) Find the equation of the tangent plane to the graph of f at the point (2,1,3).
b) Plot the graph of f and its tangent plane at (2,1,3).

Solution: Here, a=2, b=1.
a)

Clear[f,x,y,z]
f[x_,y_]=(x^2)+(y^2)

Thus, the equation the of the tangent plane is

A=(∂_x) f[x,y] /.{x->2,y->1};
 B=(∂_y) f[x,y] /.{x->2,y->1};
z=A (x-2) + B(y-1)+f[2,1];
Simplify[z]

b) Here is a plot of the graph of f:

plot1=Plot3D[{f[x,y],z},{x,-10,10},{y, -10, 10},PlotStyle-> {Blue,Green}];
plot2=ListPointPlot3D[{ {2,1,3}},PlotStyle-> {Red,PointSize[Large]}  ];
Show[plot1,plot2,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

Example 14.26.  Let f(x, y) =(x^2) y-6 x (y^2)+3y.  Find the points where  the tangent plane to the graph of f is parallel to the xy-plane.

Solution:  For the tangent plane to be parallel to the xy-plane, we must have (Subscript[f, x])=0 and (Subscript[f, y])=0.

Clear[f,x,y ]
f[x_,y_]= (x^(2 ))y-6 x (y^2)+3y

A tangent plane is parallel to the xy-plane at

Solve[{ D[f[x,y],x] ==0,  D[f[x,y],y] ==0}]

Rotate the following graph to see the points of tangencies.

Plot3D[{f[x,y],f[-1,-1/3],f[1,1/3]},{x,-1,1},{y,-1,1},PlotStyle-> {LightBlue,Green,Red},PlotRange->All,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

Exercises

1. Let  f(x,y)=(x^(3 ))y +x (y^2)-3x +4.
a) Find a set of parametric equations of the normal line and an equation of the tangent plane to the surface at the point (1, 2).
b) Graph the surface, the normal line, and the tangent plane found in a).

2. Let f(x, y) =  (x^2)+(y^2).
a. Find the equation of the tangent plane to the graph of f at the point (2,1,5).
b. Plot the graph of f and its tangent plane at (2,1,5).

3.  Let f(x, y) =  (e^(-y/x)).
a. Find the equation of the tangent plane to the graph of f at the point (1,0,1).
b. Plot the graph of f and its tangent plane at (1,0,1).

4.  Let f(x, y) =cos (x y).  Find the points where  the tangent plane to the graph of f is parallel to the xy-plane.

14.5  Gradient and Directional Derivatives

Students should read Section 14.5 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Recall that the notation for a vector such as u=2i +5j -6k in Mathematica is {2,5,-6}.  The command for the dot product of two vectors u and v is obtained by typing u.v.

The gradient of f, denoted by  ∇f, at (a,b) can be obtained by evaluating ∇f(a,b)=⟨(∂_x)f(a,b), (∂_y)f(a,b)⟩.

The directional derivative of f at (a,b) in the direction of a unit vector u is given by (Subscript[D, u])f=∇f(a,b)·u.

Example 14.27.  Find the gradient and directional derivative of f(x,y)=(x^2)sin 2y at the point  (1,(π/2),0)   in the direction of v=⟨(3/5),-(4/5)⟩.

Solution:

Clear[f,v]
f[x_,y_]:=(x^2) *Sin[2y]
v={(3/5),((-4)/5)}

The gradient of f at (1,(π/2))  is

∇f ={(∂_x) f[x,y],(∂_y) f[x,y]}/.{x->1,y->(π/2)}

Since v is a unit vector, the directional derivative is given by

direcderiv=v.∇f

Example 14.28.  Find the gradient and directional derivative of  f(x,y,z)=x y+y z+x z at the point  (1, 1, 1)  in the direction of v=2i+j-k.

Solution:

Clear[x,y,z]
w=x*y+y*z+x*z
v={2,1,-1}

We normalize v:

unitvector=v/Norm[v]

The gradient of w=f(x,y,z) at (1,1,1) is

∇w={D[w,x],D[w,y],D[w,z]}/.{x->1,y->1,z->1}

Hence, the directional derivative is given by

direcderiv=unitvector.∇w

Example 14.29.  Plot the gradient vector field and the level curves of the function f(x,y)=(x^2)sin 2y.

Solution:

Clear[f,fx,fy,x,y]
f[x_,y_]= (x^2)-3x y+y-(y^2)
fx=D[f[x,y],x]
fy=D[f[x,y],y]

Thus, the gradient vector field is ∇f(x,y)=⟨ 2x-3y,1-3x-2y⟩.  To plot this vector field, we need to download the package VectorFieldPlots, which is done by evaluating

Needs["VectorFieldPlots`"]

Here is a plot of some level curves and the gradient field.

Clear[plot1,plot2]
plot1=ContourPlot[f[x,y],{x, -5,5},{y,-4,4},Axes->True,Frame->False,Contours->15, ColorFunction->Hue] ;
plot2=VectorFieldPlot[{fx,fy},{x, -5,5},{y,-4,4},Axes->True,Frame->False];
Show[plot1, plot2,ImageSize->{250}]

Example 14.30.  Let the temperature T at a point (x,y) on a metal plate be given by  T(x,y)=(x/((x^2)+(y^2))). 
a) Plot the graph of the temperature. 
b) Find the rate of change of temperature at (3,4), in the direction of v=i-2j.
c) Find the unit vector in the direction of which the temperature increases most rapidly at (3,4).
d) Find the maximum rate of increase in the temperature at (3,4).

Solution: 
a)  Here is the graph of T.

T[x_,y_]=(x/((x^2)+(y^2)))

graphofT=Plot3D[T[x,y],{x,-5,5},{y,-5,5},BoxRatios->{1,1,1},ImageSize->Small]

b) Let u=(v/(‖v‖)).  Then u is a unit vector and the rate of change in temperature at (3,4) in the direction of v is given by (Subscript[D, u])T(3,4)=∇f(3,4)·u.

∇T={D[T[x,y],x],D[T[x,y],y]}
v={1,-2}
u=(v/(Sqrt[v.v]))
u.∇T/.{x->3,y->4}//N

Thus, the rate of change at (3,4) in the direction v is 0.0393548.  NOTE: The command //N in the last line of the previous input converts the output to decimal form.

c) The unit vector in the direction of which the temperature increases most rapidly at (3,4) is given by

((∇T)/(Norm[∇T]))/.{x->3,y->4}

d) The maximum rate of increase in the temperature at (3,4) is the norm of the gradient at this point. This can be obtained by:

Norm[∇T]/.{x->3,y->4}

Exercises

1. Find the gradient and directional derivative of  f(x,y)=(sin^(-1))(x y) at the point  (1,1,(π/2))   in the direction of v=⟨1,-1⟩.

2. Let  T(x,y)=(e^(x y -(y^2))).
a. Find ∇ T(x, y).
b. Find the directional derivative of  T(x,y) at the point (3,5) in the dierection of u = 1/2 i  + (Sqrt[3])/2j.
c. Find the direction of greatest increase in T from the point (3,5).

3.  Plot the gradient vector field and the level curves of the function a f( x , y ) = cos x ((sin )^2)y.

4.  Find the gradient and directional derivative of  f(x,y,z) =  x y (e^( y z))+ sin (x z) at the point (1,1,0)  in the direction of v=i-j-k.

14.6  The Chain Rule

Students should read Section 14.6 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Example 14.31.  Let x=(t^2)+s ,y=t+(s^(2 )) and  z=x sin y.  Find the first partial derivatives of z with respect to s and t.

Solution:

Clear[x,y,z,s,t]
x=(t^2)+s
y=t+(s^2)
z=x Sin[y]

D[z,s]

D[z,t]

Example 14.32.  Find the partial derivatives of z with respect to x and y assuming that the equation (x^2)z-y (z^2)=x y defines z as a function of x and y.

Solution:

Clear[x,y,z,r,t,s]
eq= (x^2)z[x,y] -y ((z[x,y])^2)==x y

Solve[D[eq,x],D[z[x,y],x]]
Solve[D[eq,y],D[z[x,y],y]]

Example 14.33.  Let f(x,y,z)=F(r), where r=(Sqrt[(x^2)+(y^2)+(z^2)]) and F is a twice differentiable function of one variable. 
a) Show that ∇f=F'(r) (1/r)(x i +y j+z k).
b) Find the Laplacian of f.

Solution: 
a)

Clear[x,y,z,r,f,F]
f[x_,y_,z_]=F[r]

r=(Sqrt[(x^2)+(y^2)+(z^2)])

Here is the gradient of f:

gradf={D[f[x,y,z],x],D[f[x,y,z],y],D[f[x,y,z],z]}

With r=(Sqrt[(x^2)+(y^2)+(z^2)]), the preceding output becomes

∇f(x,y,z)=⟨  ((x F'(r))/r), ((y F'(r))/r),((z F'(r))/r)⟩=F'(r)(1/r)⟨ x ,y ,z ⟩

which proves part a).

b) Recall that the Laplacian of f, denoted by Δ f, is defined by Δ f=(Subscript[f, xx])+(Subscript[f, yy])+(Subscript[f, zz]).

D[f[x,y,z],{x,2}]+D[f[x,y,z],{y,2}]+D[f[x,y,z],{z,2}]

We simplify this to get

Simplify[%]

which is the same as (2/r)F'[r]+F''[r].

Exercises

1. Let x=(u^2)+sin v, y=u (e^(v/u)), and z=(y^3)ln x  .  Find the first partial derivatives of z with respect to u and v.

2. Find the partial derivatives of z with respect to x and y assuming that the equation (x^2)z-y (z^2)=x y defines z as a function of x and y.

3. Find an equation of the tangent plane to the surface x z +2 (x^2)y +(y^2)(z^3)=11 at (2,1,1).

14.7  Optimization

Students should read Section 14.7 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Second Derivative Test: Suppose (Subscript[f, x])(a,b)=0 and (Subscript[f, y])(a,b)=0.  Define

D(x,y)=(Subscript[f, x x])(Subscript[f, y y ])-((((Subscript[f, x y])))^2)

The function D is called the discriminant function.

i)   If D(a,,b)>0 and (Subscript[f, x x])(a,b)>0, then f(a,b) is a local minimum value.
ii)   If D(a,,b)>0 and (Subscript[f, x x])(a,b)<0, then f(a,b) is a local maximum value.
iii)  If D(a,,b)<0, then (a,b,f(a,b)) is a saddle point on the graph of f.
iv)  If D(a,b)=0, then no conclusion can be drawn about the the point (a,b).

Example 14.34.  Let f(x,y)=(x^4)-4 x y+2(y^2). 
a)  Find all critical points of f.
b)  Use the second derivative test to classify the critical points as local minimum, local maximum, saddle point, or neither.

Solution: Since D is used in Mathematica as the command for derivative, we will use disc for the discriminant function D.

Clear[f,x,y]
f[x_,y_]=(x^4)-4x y +2(y^2)

a) The critical points are given by

cp=Solve[{D[f[x,y],x]==0,D[f[x,y],y]==0}]

b)

Clear[fxx,disc]
fxx[x_,y_]=D[f[x,y],{x,2}]
disc[x_,y_]=D[f[x,y],{x,2}]*D[f[x,y],{y,2}]-(((D[D[f[x,y],x],y]))^2)

TableForm[Table[{ cp[[k,2,2]],   cp[[k,1,2]] ,  disc[cp[[k,2,2]],cp[[k,1,2]]], fxx[cp[[k,2,2]],cp[[k,1,2]]],f[cp[[k,2,2]],cp[[k,1,2]]]},{k,1,Length[cp]}],
TableHeadings->{{},{"x ", "y  ",  "  D(x,y) ","  (Subscript[f, xx])  ","f(x,y)"}}]

By the second derivative test, we conclude that f has a local minimum value of -1 at (-1,-1) and (1,1), and a saddle point at (0,0).

Here is the graph of f and the relevant points.

Clear[plot1,plot2]
plot1=Plot3D[f[x,y],{x,-2,2},{y, -2,2},PlotStyle-> LightBlue,PlotRange-> 10];
plot2=Graphics3D[{PointSize[Large], Red, Point[Table[{ cp[[k,2,2]],   cp[[k,1,2]] ,  f[cp[[k,2,2]],cp[[k,1,2]]]},{k,1,Length[cp]}]]}, PlotRange->10];
Show[plot1,plot2,ImageSize->{250}]

Example 14.35.  Let f(x,y)=(x^3)+(y^4)-6x-2(y^2). 
a)  Find all critical points of f.
b)  Use the second derivative test to classify the critical points as local minimum, local maximum, saddle point, or neither.

Solution: Again, we will use disc to denote the discriminant function D since the letter D is used in Mathematica for the derivative command.

Clear[f,x,y]
f[x_,y_]=(x^3)+(y^4)-6 x -2(y^2)

a) The critical points are given by

cp=Solve[{D[f[x,y],x]==0,D[f[x,y],y]==0}]

b)

Clear[fxx,disc]
fxx[x_,y_]=D[f[x,y],{x,2}]
disc[x_,y_]=D[f[x,y],{x,2}]*D[f[x,y],{y,2}]-(((D[D[f[x,y],x],y]))^2)

TableForm[Table[{ cp[[k,2,2]],   cp[[k,1,2]] ,  disc[cp[[k,2,2]],cp[[k,1,2]]], fxx[cp[[k,2,2]],cp[[k,1,2]]],f[cp[[k,2,2]],cp[[k,1,2]]]},{k,1,Length[cp]}],
TableHeadings->{{},{"x ", "y  ",  "  D(x,y) ","  (Subscript[f, xx])  ","f(x,y)"}}]

By the second derivative test we conclude that f has local maximum value of 4(Sqrt[2]) at (-(Sqrt[2]),0), local minimum value of -1-4(Sqrt[2]) at ((Sqrt[2]),-1) and ((Sqrt[2]),1), and saddle points at (-(Sqrt[2]),-1), ((Sqrt[2]),0), and  (-(Sqrt[2]),1).

Here is the graph of f and the relevant points.

Clear[plot1,plot2]
plot1=Plot3D[f[x,y],{x,-2.5,2.5},{y, -2.5, 2.5},PlotStyle-> LightBlue,PlotRange-> 10];
plot2=Graphics3D[{PointSize[Large], Red, Point[Table[{ cp[[k,2,2]],   cp[[k,1,2]] ,  f[cp[[k,2,2]],cp[[k,1,2]]]},{k,1,Length[cp]}]]}, PlotRange->10];
Show[plot1,plot2,ImageSize->{250}]

Example 14.36.  Let f(x,y)=2(x^2)-3x y-x+y+(y^2) and let R be the rectangle in the xy-plane whose vertices are at (0,0), (2,0), (2,2), and  (0,2). 
a) Find all relative extreme values of f inside R.
b) Find the maximum and minimum values of f on R.

Solution:

Clear[f,x,y,disc]
f[x_,y_]=2(x^2)-3x*y-x+y+(y^2)+5

Solve[{(∂_x) f[x,y]==0,(∂_y) f[x,y]==0},{x,y}]

disc[x_,y_]= (∂_(x,x))f[x,y]*(∂_(y,y))f[x,y]-((((∂_(x,y))f[x,y]))^2)

(∂_(x,x))f[x,y]/.{x->1,y->1}
disc[x,y]/.{x->1,y->1}

Thus, (1,1) is the local minimum point of f inside R and its local minimum value is f(1,1)=5. Next, we find the extreme values of f on the boundary of the rectangle. This is done by considering f as a function of one variable corresponding to each side of R.  Let (Subscript[f, 1])=f(x,0), (Subscript[f, 2])=f(x,2), for x between 0 and 2, and (Subscript[f, 3])=f(0,y) and (Subscript[f, 4])=f(2,y), for y between 0 and 2.   We now proceed as follows:

Clear[f1,f2,f3,f4]
f1=f[x,0]
f2=f[x,2]
f3=f[0,y]
f4=f[2,y]

Solve[D[f1,x]==0 ]

Solve[D[f2,x]==0 ]

Solve[D[f3,y]==0 ]

Solve[D[f4,y]==0 ]

Thus, points on the boundary of R that are critical points of f are ((1/4),  0) and ((7/4), 2).  Observe that the points (0, -1/2) and (2,(5/2)) are outside the rectangle R. The four vertices of R at (0,0), (2,0), (0,2) and (2,2) are also critical points.  Can you explain why?  We now evaluate f at each of these points and at (1,1) (the relative minimum point found earlier) using the substitution command and compare the results.

f[x,y]/.{{x->(1/4),y->0},  {x->(7/4),y->2},
		{x->0,y->0},{x->2,y->0},{x->0,y->2},{x->2,y->2},{x->1,y->1}}

Thus, the minimum value of f is 39/8, which occurs at (1/4,0) and also at (7/4,2).  The maximum value of f is 6, which is attained at (2,0) and also at (0,2).  Here is the graph of f over the rectangle R.

Clear[plot1,plot2,plot3]
plot1=Plot3D[{f[x,y],0},{x,0,2},{y,0,2},PlotStyle->{Green,Blue},PlotRange->All];
plot2=Graphics3D[{PointSize[Large], Red , Point[{ {1/4,0,f[1/4,0]},{7/4,2,f[7/4,2]} }]}, PlotRange->All ];
plot3=Graphics3D[{PointSize[Large],Black , Point[{  {2,0,f[2,0]},{0,2,f[2,0]}}]}, PlotRange->All ];
Show[plot1,plot2 ,plot3,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

Exercises

1. Let f(x,y)=(x^4)-4 x y+2(y^2). 
a. Find all critical points of f.
b. Use the second derivative test to classify the critical points as local minimum, local maximum, saddle point, or neither.
c. Plot the graph of f and the local extreme points and saddle points, if any.

2. Let f(x,y)= (x+y)ln((x^2)+(y^2)), for (x,y)≠(0,0). 
a. Find all critical points of f.
b. Use the second derivative test to classify the critical points as local minimum, local maximum, saddle point, or neither.
c. Plot the graph of f and the local extreme points and saddle points, if any.

3. Let f(x,y)=2(x^2)-3x y-x+y+(y^2) and let R be the rectangle in the xy-plane whose vertices are at (0,0), (2,0), (2,2), and (0,2). 
a. Find all relative extreme values of f inside R.
b. Find the maximum and minimum values of f on  R.
c. Plot the graph of f and the local extreme points and saddle points, if any.

14.8  Lagrange Multipliers

Students should read Section 14.8 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Example 14.37.   Let f(x,y)=x y and g(x,y)=(x^2)+(y^2)-4.
a) Plot the level curves of f and g as well as their gradient vectors. 
b) Find the maximum and minimum values of f subject to the constraint g(x,y)=0.

Solution: 
a) We will define f and g and compute their gradients. Recall that we need to evaluate the command Needs["`VectorFieldPlots`"] before we plot the gradient fields.

Clear[f,g,fx,fy,gx,gy,x,y]
f[x_,y_]=2x+3y
g[x_,y_]=(x^2)+(y^2)-4
fx=D[f[x,y],x]
fy=D[f[x,y],y]
gx=D[g[x,y],x]
gy=D[g[x,y],y]

Needs["VectorFieldPlots`"]

Clear[plot1,plot2,plot3,plot4]
plot1=ContourPlot[(x^2)+(y^2)-4,{x,-2,2},{y,-2,2},Frame->False,Axes->True, ContourShading->False,PlotRange->All];
plot2= ContourPlot[2x+3y,{x,-2,2},{y,-2,2},Frame->False,Axes->True, ContourShading->False,PlotRange->All];
plot3=VectorFieldPlot[{fx,fy},{x,-2,2},{y,-2,2},Axes->True,Frame->False, ColorFunction->Hue];
plot4=VectorFieldPlot[{gx,gy},{x,-2,2},{y,-2,2},Axes->True,Frame->False, ColorFunction->Hue];
 Show[plot1, plot2,plot3,plot4,ImageSize->{250}]

b) Let us use l for λ. To solve ∇f=l ∇g we compute

Solve[{fx==  l gx,fy==l gy, g[x,y]==0} ]

Thus, (-(4/(Sqrt[13])),-(6/(Sqrt[13]))) and ((4/(Sqrt[13])),(6/(Sqrt[13]))) are the critical points. We evaluate f at these points to determine the absolute maximum and the absolute minimum of f on the graph of g(x,y)=0.

f[ -(4/(Sqrt[13])),-(6/(Sqrt[13]))]
f[(4/(Sqrt[13])),(6/(Sqrt[13]))]

Hence, f attains its absolute minimum value of -2(Sqrt[13]) at (-(4/(Sqrt[13])),-(6/(Sqrt[13]))) and absolute maximum value of -2(Sqrt[13]) at ((4/(Sqrt[13])),(6/(Sqrt[13]))).

Here is a combined plot of the gradients of f (in black) and g (in red) at the critical points.

Clear[plot1,plot2,plot3, plot4,plot5,plot6]
plot1= ContourPlot[g[x,y],{x,-3,3}, {y, -3,3}, Contours->{0},Frame->False,Axes->True, ContourShading->False];
plot2=ListPlot[{{ -(4/(Sqrt[13])),-(6/(Sqrt[13]))}, {(4/(Sqrt[13])),(6/(Sqrt[13]))}}];

plot3=Graphics[Arrow[{{ -(4/(Sqrt[13])),-(6/(Sqrt[13]))},{ -(4/(Sqrt[13])),-(6/(Sqrt[13]))} +{fx,fy}/.{x->  ((-4)/(Sqrt[13])),y->((-6)/(Sqrt[13]))}  }]];

plot4=Graphics[Arrow[{ {(4/(Sqrt[13])),(6/(Sqrt[13]))} ,{(4/(Sqrt[13])),(6/(Sqrt[13]))}+{fx,fy}/.{x->  (4/(Sqrt[13])),y->(6/(Sqrt[13]))}  }]];

plot5=Graphics[{Red,Arrow[{{ -(4/(Sqrt[13])),-(6/(Sqrt[13]))} ,{ -(4/(Sqrt[13])),-(6/(Sqrt[13]))}+ {gx,gy}/.{x->  ((-4)/(Sqrt[13])),y->((-6)/(Sqrt[13]))} }]}];

plot6=Graphics[{Red,Arrow[{ {(4/(Sqrt[13])),(6/(Sqrt[13]))} ,{(4/(Sqrt[13])),(6/(Sqrt[13]))}+{gx,gy}/.{x->  (4/(Sqrt[13])),y->(6/(Sqrt[13]))} }]}];

Show[plot1,plot2,plot3,plot4, plot5,plot6,PlotRange->All,AspectRatio->Automatic,ImageSize->{250}]

Exercises

1.  Let f(x,y)=4(x^2)+9(y^2) and g(x,y)=x y-4.
a. Plot the level curves of f and g as well as their gradient vectors. 
b. Find the maximum and minimum values of f subject to g(x,y)=0.

2.  Find the maximum and minimum values of f(x,y,z)=(x^3)-3(y^2)+4z subject to the constraint g(x,y,z)=x +y z-4=0.

3.  Find the maximum area of a rectangle that can be inscribed in the ellipse ((x^2)/(a^2))+((y^2)/(b^2))=1.

4.  Find the maximum volume of a box that can be inscribed in the sphere (x^2)+(y^2)+(z^2)=4.

Chapter 15	Multiple Integration

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

15.1  Double Integral over a Rectangle

Students should read Section 15.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Integration can be generalized to functions of two or more variables.  As the integral of a single-variable function defines area of a plane region under the curve, it is natural to consider a double integral of a two-variable function that defines volume of a solid under a surface.  This definition can be made precise in terms of double Riemann sums where rectangular columns (as opposed to rectangles) are used as building blocks to approximate volume (as opposed to area).  The exact volume is then obtained as a limit where the number of columns increases without bound.

15.1.1  Double Integrals and Riemann Sums

Let f(x,y) be a function of two variables defined on a rectangular domain R=[a,b]×[c,d] in (ℝ^2).  Let P={a=(Subscript[x, 0])<(Subscript[x, 1])<...<(Subscript[x, m])=b,c=(Subscript[y, 0])<(Subscript[y, 1])<...<(Subscript[y, n])=d} be an arbitrary partition of R into a grid of m·n rectangles, where m and n are integers.  For each sub-rectangle (Subscript[R, ij])=[(Subscript[x, i-1]),(Subscript[x, i])]×[(Subscript[y, j-1]),(Subscript[y, j])] denote by (Subscript[ΔA, ij]) its area and choose an arbitrary base point ((Subscript[x, ij]),(Subscript[y, ij]))∈(Subscript[R, ij]), where (Subscript[x, ij])∈[(Subscript[x, i-1]),(Subscript[x, i])] and (Subscript[y, ij])∈[(Subscript[y, j-1]),(Subscript[y, j])].  The product f((Subscript[x, ij]),(Subscript[y, ij]))(Subscript[ΔA, ij]) represents the volume of the ij-rectangular column situated between the surface and the xy-plane.  We then define the double Riemann sum (Subscript[S, p]) of f(x,y) on R with respect to P to be the total volume of all these columns:

(Subscript[S, P])=(∑_(i=1)^m)(∑_(j=1)^n)f((Subscript[x, ij]),(Subscript[y, ij]))(Subscript[ΔA, ij])

Define ||P|| to be the maximum dimension of all the sub-rectangles.  The double integral of f(x,y) on the rectangle R is then defined as the limit of (Subscript[S, P]) as ||P||->0:

∫(∫_R)f(x,y)dA=(lim_(||P||->0))(∑_(i=1)^m)(∑_(j=1)^n)f((Subscript[x, ij]),(Subscript[y, ij]))(Subscript[ΔA, ij])

If the limit exists regardless of the choice of partition and base points, then the double integral is said to exist.  Otherwise, the double integral does not exist.

MIDPOINT RULE (Uniform Partitions): Let us consider uniform partitions P, where the points {(Subscript[x, i])} and {(Subscript[y, j])} are evenly spaced, that is,  (Subscript[x, i])=a+i Δ x, (Subscript[y, j])=b+j Δ y for i=0,1,...,m and j=0,1,...,n, and with Δx=(b-a)/m and Δy=(d-c)/n.  Then the corresponding double Riemann sum is

(Subscript[S, m,n])=(∑_(i=1)^m)(∑_(j=1)^n)f((Subscript[x, ij]),(Subscript[y, ij]))Δ x Δ y

Here is a subroutine called MDOUBLERSUM that calculates the double Riemann sum (Subscript[S, m,n]) of f(x,y) over a rectangle R for uniform partitions using the center midpoint of each sub-rectangle as base point, that is, (Subscript[x, ij])=((Subscript[x, i-1])+(Subscript[x, i]))/2=a+(i-1/2)Δ x and (Subscript[y, ij])=((Subscript[y, j-1])+(Subscript[y, j]))/2=c+(j-1/2)Δ y.

Clear[f]
MDOUBLERSUM[a_,b_,c_,d_,m_,n_]:=Sum[f[a+(i-1/2)*(b-a)/m,c+(j-1/2)*(d-c)/n]*(b-a)/m*(d-c)/n,{i,1,m},{j,1,n}]

Example 15.1.  Approximate the volume of the solid bounded below the surface f(x)=(x^2)+(y^2) and above the rectangle R=[-1,1]×[-1,1] on the xy-plane using a uniform partition with m=10 and n=10 and center midpoints as base points.  Then experiment with larger values of m and n and conjecture an answer for the exact volume.

Solution: We calculate the approximate volume for m=10 and n=10 using the subroutine MDOUBLERSUM:

f[x_,y_]=x^2+y^2;
MDOUBLERSUM[-1,1,-1,1,10,10]

N[%]

Table[MDOUBLERSUM[-1,1,-1,1,10*k,10*k],{k,1,10}]

N[%]

It appears that the exact volume is 8/3.  To prove this, we evaluate the double Riemann sum (Subscript[S, m,n]) in the limit as m,n->∞:

Clear[S,m,n];
S[m_,n_]=Simplify[MDOUBLERSUM[-1,1,-1,1,m,n]]

Limit[Limit[S[m,n],m->Infinity],n->Infinity]

To see this limiting process visually, evaluate the following subroutine, called DOUBLEMIDPT, which plots the surface of the function corresponding to the double integral along with the rectangular columns defined by the double Riemann sum considered in the previous subroutine MDOUBLERSUM.

Clear[f];
DOUBLEMIDPT[f_,{a_,b_,m_},{c_,d_,n_}]:=Module[
{dx,dy, i,j, xstar,ystar,mrect, plot},
dx=N[(b-a)/m];
xstar=Table[a+i* dx, {i,0,m}]; 
dy=N[(d-c)/n];
ystar=Table[c+j* dy, {j,0,n}]; 
mcolumn=Table[Cuboid[{xstar[[i]],ystar[[j]],0}, {xstar[[i+1]],ystar[[j+1]],f[(xstar[[i]]+xstar[[i+1]])/2,(ystar[[j]]+ystar[[j+1]])/2]}],{i,1,m},{j,1,n}];
plot=Plot3D[f[x,y],{x,a,b},{y,c,d},Filling ->Bottom];
Show[plot,Graphics3D[mcolumn],ImageSize->{300}]]

f[x_,y_]:=(x^2)+y^2;
DOUBLEMIDPT[f,{-1,1,10},{-1,1,10}]

Here is an animation that demonstrates how the volume of the rectangular columns approach that of the solid in the limit as m,n->∞:

Important Note: If you are reading the printed version of this publication, then you will not be able to view any of the animations generated from the Animate command in this chapter.  If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then evaluate each Animate command to view the corresponding animation.  Just click on the arrow button to start the animation.  To control the animation just click at various points on the sliding bar or else manually drag the bar.

Animate[DOUBLEMIDPT[f,{-1,1,a},{-1,1,a}] ,{a,5,50,5 }]

15.1.2  Double Integrals and Iterated Integrals in Mathematica

The Mathematica command for evaluating double integrals is the same as that for evaluating integrals of a single-variable function, except that two limits of integration must be specified, one for each independent variable.  Thus:

Integrate[f[x,y],{x,a,c},{y,c,d}] analytically evaluates the double integral ∫(∫_R)f(x,y)dA over the rectangle R=[a,b]×[c,d].

NIntegrate[f[x,y],{x,a,c},{y,c,d}] numerically evaluates the double integral ∫(∫_R)f(x,y)dA over the rectangle R=[a,b]×[c,d].

Iterated Integrals: In practice, one does not actually use the limit definition in terms of Riemann sums to evaluate double integrals, but instead apply Fubini's Theorem to easily compute them in terms of iterated integrals:

Fubini's Theorem: (Rectangular Domains) If R={(x,y):a≤x≤b,c≤y≤d}, then

∫(∫_R)f(x,y)dA=(∫_a^b)(∫_c^d)f(x,y)dydx=(∫_c^d)(∫_a^b)f(x,y)dxdy

Thus, Mathematica will naturally apply Fubini's Theorem whenever possible to analytically determine the answer.  Depending on the form of the double integral, Mathematica may resort to more sophisticated integration techniques, such as contour integration, which are beyond the scope of this text.

Example 15.2.  Calculate the volume of the solid bounded below by the surface f(x)=(x^2)+(y^2) and above the rectangle R=[-1,1]×[-1,1].

Solution: The volume of the solid is given by the double integral ∫(∫_R)f(x,y)dA.  To evaluate it, we use the Integrate command:

f[x_,y_]:=x^2+y^2;
Integrate[f[x,y],{x,-1,1},{y,-1,1}]

This confirms the conjecture that we made in the previous example for the exact volume.

NOTE: Observe that we obtain the same answer by explicitly computing this double integral as an integrated integral as follows.  Moreover, for rectangular domains, the order of integration does not matter.

Integrate[Integrate[f[x,y],{x,-1,1}],{y,-1,1}]
Integrate[Integrate[f[x,y],{y,-1,1}],{x,-1,1}]

Example 15.3.  Compute the double integral ∫(∫_R)x (e^(-(y^2)))dA on the rectangle R=[0,1]×[0,1].

Solution: Observe that the Integrate command here gives us an answer in terms of the non-elementary error function Erf:

Integrate[x*E^(-y^2),{x,0,1},{y,0,1}]

This is because the function f(x,y)=x (e^(-(y^2))) has no elementary anti-derivative with respect to y due to the Gaussian factor (e^(-(y^2))) (bell curve).  Thus, we instead use the NIntegrate Command to numerically approximate the double integral:

NIntegrate[x*E^(-y^2),{x,0,1},{y,0,1}]

Exercises

1. Consider the function f(x,y)=16-(x^2)-(y^2) defined over the rectangle R=[0,2]×[-1,3].
a. Use the subroutine MDOUBLERSUM to compute the double Riemann sum (Subscript[S, m,n]) of f(x,y) over R for m=2 and n=2.
b. Repeat part a) by generating a table of double Riemann sums for m=10k and n=10k where k=1,2,...,10.  Make a conjecture for the exact value of ∫(∫_R)f(x,y)dA. 
c. Find a formula for (Subscript[S, m,n]) in terms of m and n.  Verify your conjecture in part b) by evaluating (lim_(m,n->∞))(Subscript[S, m,n]).
d. Directly compute ∫(∫_R)f(x,y)dA using the Integrate command.

2. Repeat Exercise 1 but with f(x,y)=(1+x)(1+y)(1+x y) defined over the rectangle R=[0,1]×[0,1].

3. Evaluate the double integral ∫∫(Sqrt[(x^4)+(y^4)])dA over the rectangle R=[-2,1]×[-1,2] using both the Integrate and NIntegrate commands.  How do the two answers compare?

4. Calculate the volume of the solid lying under the surface z=(e^(-y))(x+(y^2)) and over the rectangle R=[0,2]×[0,3].  Then make a plot of this solid.

5. Repeat Exercise 4 but with z=sin((x^2)+(y^2)) and rectangle R=[-(Sqrt[π]),(Sqrt[π])]×[-(Sqrt[π]),(Sqrt[π])].

6. Evaluate the double integral ∫(∫_R)f(x,y)dA where f(x,y)=x y cos((x^2)+(y^2)) and R=[-π,π]×[-π,π].  Does your answer make sense?  Make a plot of the solid corresponding to this double integral to intuitively explain your answer.  HINT: Consider symmetry.

7. Find the volume of solid bounded between the two hyperbolic paraboloids (saddles) z=1+(x^2)-(y^2) and z=3-(x^2)+(y^2) over the rectangle R = [-1,1]×[-1,1].

8. Find the volume of the solid bounded by the planes z=2x, z=-3x+2, y=0, y=1, and z=0.

15.2  Double Integral over More General Regions

Students should read Section 15.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

For domains of integration that are non-rectangular but still simple, that is, bounded between two curves, Fubini's Theorem continues to hold.  There are two types to consider:

Fubini's Theorem: (Simple Domains)

Type I (Vertically Simple): If D={(x,y):a≤x≤b,α(x)≤y≤β(x)}, then

∫(∫_D)f(x,y)dA=(∫_a^b)(∫_(α(x))^(β(x)))f(x,y)dydx

The corresponding Mathematica command is Integrate[f[x,y],{x,a,b},{y,α[x],β[x]}].

Type II (Horizontally Simple): If D={(x,y):c≤y≤d,α(y)≤x≤β(y)}, then

∫(∫_D)f(x,y)dA=(∫_c^d)(∫_(α(y))^(β(y)))f(x,y)dxdy

The corresponding Mathematica command is Integrate[f[x,y],{y,c,d},{x,α[y],β[y]}].

Warning: Be careful not to reverse the order of integration prescribed for either type.  For example, evaluating the command Integrate[f[x,y],{y,α[x],β[x]},{x,a,b}] for Type I (x and y are reversed) will lead to incorrect results.

Example 15.4.  Calculate the volume of the solid bounded below by the surface f(x,y)=1-(x^2)+(y^2) and above the domain D bounded by x=0, x=1, y=x, and y=1+(x^2).

Solution: We observe that x=0 and x=1 represent the left and right boundaries, respectively, of D.  Therefore, we plot the graphs of the other two equations along the x-interval [0,2] to visualize D (shaded in the following plot):

Clear[x,y]
plot1=Plot[{x,1+x^2},{x,0,1},Filling->{1->{2}},ImageSize->250]

Here is a plot of the corresponding solid situated over D:

f[x_,y_]=1-x^2+y^2;
plot3=Plot3D[f[x,y],{x,0,1},{y,x,1+x^2},Filling->Bottom,ViewPoint->{1,1,1},PlotRange->{0,4},ImageSize->{250}]

To compute the volume of this solid given by ∫(∫_D)f(x,y)dA, we describe D as a vertically simple domain where 0≤x≤1 and x≤y≤1+(x^2) and apply Fubini's Theorem to evaluate the corresponding iterated integral (∫_0^1)(∫_x^(1+(x^2)))f(x,y)dydx (remember to use the correct order of integration):

Integrate[f[x,y],{x,0,1},{y,x,1+(x^2)}]

Example 15.5.  Evaluate the double integral ∫(∫_D)sin ((y^2))dA, where D is the domain bounded by x=0, y=2, and y=x.

Solution: We first plot the graphs of  x=0, y=2, and y=x to visualize the domain D:

plot1=ContourPlot[{x==0,y==2,y==x},{x,-0.5,2.5},{y,-0.5,2.5},ImageSize->{250}]

It follows that D is the triangular region bounded by these graphs, which we shade in the following plot to make clear:

plot2=Plot[x,{x,0,2},Filling->2];
Show[plot1,plot2,ImageSize->{250}]

To compute the given double integral, we describe D as a horizontally simple domain, where 0≤y≤2 and 0≤x≤y and apply Fubini's Theorem to evaluate the corresponding iterated integral (∫_0^2)(∫_0^y)sin((y^2))dxdy (again, remember to use the correct order of integration):

Integrate[Sin[y^2],{y,0,2},{x,0,y}]

N[%]

NOTE: It is also possible to view D as a vertically simple domain, where 0≤x≤2 and x≤y≤2.  The corresponding iterated integral (∫_0^2)(∫_x^2)sin((y^2))dydx gives the same answer, as it should by Fubini's Theorem:

Integrate[Sin[y^2],{x,0,2},{y,x,2}]

Observe that it is actually impossible to evaluate this iterated integral by hand since there is no elementary formula for the anti-derivative of sin((y^2)) with respect to y.  Thus, if necessary, Mathematica automatically switches the order of integration by converting from one type to the other.

Exercises

In Exercises 1 through 4, evaluate the given iterated integrals and plot the solid corresponding to each one.
1. (∫_0^1)(∫_0^(x^2))(4-(x^2)+(y^2))dydx		2. (∫_0^4)(∫_0^(Sqrt[2-(y^2)]))(x^2) ydxdy
3. (∫_0^π)(∫_0^(sin θ))(r^2)cosθdrdθ		4. (∫_0^1)(∫_0^x)(y/(1+x y))dydx

In Exercises 5 through 8, evaluate the given double integrals and plot the solid corresponding to each one.
5. ∫(∫_D)(x+y)dA,  D={(x,y):0≤x≤3,0≤y≤(Sqrt[x])}
6. ∫(∫_D)(Sqrt[x+y])dA,  D={(x,y):0≤x≤1-(y^2),0≤y≤1}
7. ∫(∫_D)(e^(x+y))dA, where D={(x,y):(x^2)+(y^2)≤4}
8. ∫(∫_D)(y/(x+1))dA, where D is the following shaded diamond region:

In Exercises 9 through 12, calculate the volume of the given solid S:
9. S is bounded under the paraboloid z=16-(x^2)-(y^2) and above the region bounded between the line y=x and the parabola y=6-(x^2).
10. S is bounded under the right circular cone z=(Sqrt[(x^2)+(y^2)]) and above the disk (x^2)+(y^2)≤1.
11. S is bounbed between the plane z=5+2x+2y and the paraboloid z=12-(x^2)-(y^2).  HINT: Equate the two surfaces to obtain the equation of the domain.

12. S is bounded between the cylinders (x^2)+(y^2)=1 and (y^2)+(z^2)=1.

15.3  Triple Integrals

Students should read Section 15.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Once the notion of a double integral is well established, it is straightforward to generalize it to triple (and even higher-order) integrals for functions of three variables defined over a solid region in space.  Here is the definition of a triple integral in terms of triple Riemann sums for a function f(x,y,z) defined on a box region B={(x,y,z):a≤x≤b,c≤y≤d,p≤z≤q} (refer to your calculus text for details):

∫∫(∫_B)f(x,y)dV=(lim_(||P||->∞))(∑_(i=1)^m)(∑_(j=1)^n)(∑_(k=1)^p)f((Subscript[x, ijk]),(Subscript[y, ijk]))(Subscript[ΔV, ijk])

where the notation is analogous to that used for double integrals in Section 15.1 of this text.  Of course, Fubini's Theorem also generalizes to triple integrals:

Fubini's Theorem: (Box Domains) If B={(x,y,z):a≤x≤b,c≤y≤d,p≤z≤q}, then

∫∫(∫_B)f(x,y)dV=(∫_a^b)(∫_c^d)(∫_p^q)f(x,y)dzdydx

The corresponding Mathematica commands are:

Integrate[f[x,y,z],{x,a,c},{y,c,d},{z,e,f}] analytically evaluates the triple integral ∫∫(∫_B)f(x,y)dV over the box B=[a,b]×[c,d]×[e,f].

NIntegrate[f[x,y],{x,a,c},{y,c,d},{z,e,f}] numerically evaluates the triple integral ∫∫(∫_B)f(x,y)dV over the rectangle B=[a,b]×[c,d]×[e,f].

NOTE: For box domains, the order of integration does not matter so that it is possible to write five other versions of triple iterated integrals besides the one given in Fubini's Theorem.

Example 15.6.  Calculate the triple integral  ∫∫(∫_B)x y zdV over the box B=[0,1]×[2,3]×[4,5].

Solution: We use the Integrate command to calculate the given triple integral.

Integrate[x y z,{x,0,1},{y,2,3},{z,4,5}]

Volume as Triple Integral: Recall that if a solid region W is bounded between two surfaces ψ(x,y) and ϕ(x,y), where both are defined on the same domain D with ψ(x,y)≤ϕ(x,y), then its volume V can be expressed by the triple integral

V=∫∫(∫_W)1dV=∫(∫_D)(∫_(ψ(x,y))^(ϕ(x,y)))1dzdA

Example 15.7.  Calculate the volume of the solid bounded between the surfaces z=4(x^2)+4(y^2) and z=16-4(x^2)-4(y^2) on the rectangular domain [-1,1]×[-1,1].

Solution: Here is a plot of the solid:

Plot3D[{4x^2+4y^2,16-4x^2-4y^2},{x,-1,1},{y,-1,1},Filling ->{1->8,2->8},ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

The volume of the solid is given by the triple iterated integral (∫_(-1)^1)(∫_(-1)^1)(∫_(4(x^2)+4(y^2))^(16-4(x^2)-4(y^2)))1dzdydx:

Integrate[1,{x,0,1},{y,-1,0},{z,4x^2+4y^2,16-4x^2-y^2}]

Exercises

In Exercises 1 through 4, evaluate the given iterated integrals:
1. (∫_0^1)(∫_0^x)(∫_0^(y^2))(x+y+z)dzdydx		2. (∫_0^3)(∫_0^(sin y))(∫_0^(y +z))x y zdxdzdy 
3. (∫_0^π)(∫_0^θ)(∫_0^(r cos θ))r (z^2)dzdrdθ			4. (∫_0^1)(∫_(Sqrt[x])^1)(∫_0^(1-y))z log(1+x y)dzdydx

In Exercises 5 through 8, evaluate the given triple integrals:
5. ∫∫(∫_W)(x+y z)dV, where W={(x,y,z):0≤x≤1,0≤y≤(Sqrt[x]),0≤z≤(y^2)}.
6. ∫∫(∫_W)sin ydV, where W lies under the plane z=1+x+y and above the triangular region bounded by x=0, x=2, and y=3x. 
7. ∫∫(∫_W)zdV, where W is bounded by the paraboloid z=4-(x^2)-(y^2) and z=0.
8. ∫∫(∫_W)f(x,y,z)dV, where f(x,y,z)=(z^2) and W is bounded between the cone z=(Sqrt[(x^2)+(y^2)]) and z=9.

9. The triple integral (∫_0^1)(∫_(x/2)^(1-x/2))(∫_0^(2-x-2y))dzdydx represents the volume of a solid S.  Evaluate this integral.  Then make a plot of S and describe it.

10. Midpoint Rule for Triple Integrals:
a. Develop a subroutine called MTRIPLERSUM to compute the triple Riemann sum of the triple integral ∫∫(∫_B)f(x,y,z)dV over the box domain B={(x,y,z):a≤x≤b,c≤y≤d,p≤z≤q} for uniform partitions and using the center midpoint of each sub-box as base point.  HINT: Modify the subroutine MDOUBLERSUM in Section 15.1 of this text.

b. Use your subroutine MTRIPLESUM in part a) to compute the triple Riemann sum of ∫∫(∫_B)((((x^2)+(y^2)+(z^2)))^(3/2))dV over the box B={(x,y,z):0≤x≤1,0≤y≤2,0≤z≤3} by dividing B into 48 equal sub-boxes, that is, cubes having side length of 1/2.

c. Repeat part b) by dividing B into cubes having side length of 1/4 and more generally into cubes having side length of 1/(2^n) for n sufficiently large in order to obtain an approximation accurate to 2 decimal places.

d. Verify your answer in part c) using Mathematica's NIntegrate command.

15.4  Integration in Polar, Cylindrical, and Spherical Coordinates

Students should read Section 15.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

15.4.1  Double Integrals in Polar Coordinates

The following Change of Variables Formula converts a double integral in rectangular coordinates to one in polar coordinates:

Change of Variables Formula (Polar Coordinates):

I. Polar Rectangles: If R={(r,θ):(Subscript[θ, 1])≤θ≤(Subscript[θ, 2]),(Subscript[r, 1])≤r≤(Subscript[r, 2])}, then

∫(∫_R)f(x,y)dA=(∫_(Subscript[θ, 1])^(Subscript[θ, 2]))(∫_(Subscript[r, 1])^(Subscript[r, 2]))f(r cos θ,r sin θ)rdrdθ

II. Polar Regions: If D={(r,θ):(Subscript[θ, 1])≤θ≤(Subscript[θ, 2]),α(θ)≤r≤β(θ)}, then

∫(∫_D)f(x,y)dA=(∫_(Subscript[θ, 1])^(Subscript[θ, 2]))(∫_(α(θ))^(β(θ)))f(r cos θ,r sin θ)rdrdθ

Example 15.8.  Calculate the volume of the solid region bounded by the paraboloid f(x)=4-(x^2)-(y^2) and the xy-plane using polar coordinates.

Solution: We first plot the paraboloid:

f[x_,y_]=4-x^2-y^2
Plot3D[f[x,y],{x,-2,2},{y,-2,2},PlotRange->{0,4},ImageSize->{250}]

The circular domain D can be easily described in polar coordinates by the polar rectangle R={(r,θ):0≤r≤2,0≤θ≤2π}.  Thus, the volume of the solid is given by the corresponding double integral (∫_0^(2π))(∫_0^2)f(r cos θ,r sin θ)rdrdθ in polar coordinates:

Clear[r,θ];
Integrate[r*f[r*Cos[θ],r*Sin[θ]],{r,0,2},{θ,0,2Pi}]

Observe that here f(x,y) simplifies nicely in polar coordinates:

f[r*Cos[θ],r*Sin[θ]]
Simplify[%]

NOTE: Evaluating the same double integral in rectangular coordinates by hand would be quite tedious.  This is not a problem with Mathematica, however:

Integrate[f[x,y],{x,-2,2},{y,-Sqrt[4-x^2],Sqrt[4-x^2]}]

15.4.2  Triple Integrals in Cylindrical Coordinates

The following Change of Variables Formula converts a triple integral in rectangular coordinates to one in cylindrical coordinates:

Change of Variables Formula (Cylindrical Coordinates): If a solid region W is described by (Subscript[θ, 1])≤θ≤(Subscript[θ, 2]), α(θ)≤r≤β(θ), and (Subscript[z, 1])(r,θ)≤z≤(Subscript[z, 2])(r,θ), then

∫∫(∫_W)f(x,y,z)dV=(∫_(Subscript[θ, 1])^(Subscript[θ, 2]))(∫_(α(θ))^(β(θ)))(∫_((Subscript[z, 1])(r,θ))^((Subscript[z, 2])(r,θ)))f(r cos θ,r sin θ,z)rdzdrdθ

Example 15.9.  Use cylindrical coordinates to calculate the triple integral ∫∫(∫_W) zdV, where W is the solid region bounded above by the plane z=8-x-y, below by the paraboloid z=4-(x^2)-(y^2), and inside the cylinder (x^2)+(y^2)=4.

Solution: Since W lies inside the cylinder (x^2)+(y^2)=4, this implies that it has a circular base on the xy-plane given by the same equation, which can be described in polar coordinates by 0≤θ≤2π and 0≤r≤2.  Here is a plot of all three surfaces (plane, paraboloid, and cylinder):

plotplane=Plot3D[8-x-y,{x,-2,2},{y,-2,2}];
plotparaboloid=Plot3D[4-x^2-y^2,{x,-2,2},{y,-2,2}];
plotcylinder=ParametricPlot3D[{2*Cos[θ],2*Sin[θ],z},{θ,0,2π},{z,0,12}];
Show[plotplane,plotparaboloid,plotcylinder,PlotRange->All,ImageSize->{250}]

Since W is bounded in z by 4-(x^2)-(y^2)≤z≤8-x-y, or in cylindrical coordinates,  4-r cos θ-r sin θ≤z≤4-(r^2), it follows that the given triple integral transforms to

(∫_0^(2π))(∫_0^2)(∫_(4-(r^2))^(4-r cos θ-r sin θ))z rdzdrdθ

Evaluating this integral in Mathematica yields the answer

Integrate[z*r,{θ,0,2π},{r,0,2},{z,4-r*Cos[θ]-r*Sin[θ],8+r*Cos[θ]+r*Sin[θ]}]

15.4.3  Triple Integrals in Spherical Coordinates

The following Change of Variables Formula converts a triple integral in rectangular coordinates to one in spherical coordinates:

Change of Variables Formula (Spherical Coordinates): If a solid region W is described by (Subscript[θ, 1])≤θ≤(Subscript[θ, 2]), (Subscript[ϕ, 1])≤ϕ≤(Subscript[ϕ, 2]), and (Subscript[ρ, 1])(θ,ϕ)≤ρ≤(Subscript[ρ, 2])(θ,ϕ), then

∫∫(∫_W)f(x,y,z)dV=(∫_(Subscript[θ, 1])^(Subscript[θ, 2]))(∫_(Subscript[ϕ, 1])^(Subscript[ϕ, 2]))(∫_((Subscript[ρ, 1])(θ,ϕ))^((Subscript[ρ, 2])(θ,ϕ)))f(ρ cos θ sin ϕ,ρ sin θ sin ϕ,ρ cos ϕ)(ρ^2)sinϕdρdϕdθ

Example 15.10.  Use spherical coordinates to calculate the volume of the solid W lying inside the sphere (x^2)+(y^2)+(z^2)=z and above the cone z=(Sqrt[(x^2)+(y^2)]).

Solution: In spherical coordinates, the equation of the sphere is given by

((ρ )^2)=ρ  cos ϕ

or equivalently, ρ=cos ϕ.  Similarly, the equation of the cone transforms to

ρ cos ϕ=(Sqrt[(((ρ cos θ sin ϕ))^2)+(((ρ sin θ sin ϕ))^2)])=ρ  sin ϕ

It follows that cos ϕ=sin ϕ, or ϕ=π/4.  Therefore, the cone makes an angle of 45 degrees with respect to the z-axis, as shown in the following plot along with the top half of the sphere:

Clear[ρ]
plotcone=ParametricPlot3D[{ρ Cos[θ] Sin[Pi/4],ρ Sin[θ] Sin[Pi/4],ρ Cos[Pi/4]},{θ,0,2Pi},{ρ,0,Sqrt[2]/2}];
plotsphere=ParametricPlot3D[{Cos[ϕ] Cos[θ] Sin[ϕ],Cos[ϕ] Sin[θ] Sin[ϕ],Cos[ϕ] Cos[ϕ]},{θ,0,2Pi},{ϕ,0,Pi/4}];
Show[plotcone,plotsphere,PlotRange->All,ViewPoint->{1,1,1/4},ImageSize->{250}]

It is now clear that the solid W is described by 0≤θ≤2π, 0≤ϕ≤π/4, and 0≤ρ≤cos ϕ.  Thus, its volume is given by the triple integral

(∫_0^(2π))(∫_0^(π/4))(∫_0^(cos ϕ))(ρ^2)sin ϕdρdϕdθ

which in Mathematica evaluates to

Integrate[ρ^2*Sin[ϕ],{θ,0,2Pi},{ϕ,0,Pi/4},{ρ,0,Cos[ϕ]}]

Exercises

In Exercises 1 through 4, evaluate the given double integral by converting to polar coordinates:
1. (∫_(-1)^1)(∫_(-(Sqrt[1-(x^2)]))^(Sqrt[1-(x^2)]))(1-(x^2)-(y^2))dydx			2. (∫_0^2)(∫_0^(Sqrt[4-(x^2)]))(e^(-((x^2)+(y^2))))dydx
3. ∫(∫_D)x log ydA, where D is the annulus (donut-shaped region) with inner radius 1 and outer radius 3.
4. ∫(∫_D)arctan(y/x)dA, where D is the region inside the cardioid r=1+cos t.

5. Use polar coordinates to calculate the volume of the solid that lies below the paraboloid z=(x^2)+(y^2) and inside the cylinder (x^2)+(y^2)=2y.

6. Evaluate the triple integral (∫_0^2)(∫_0^(Sqrt[4-(x^2)]))(∫_0^(Sqrt[4-(x^2)-(y^2)]))((x^2)+(y^2))dzdydx by converting to cylindrical coordinates.

7. Use cylindrical coordinates to calculate the triple integral ∫∫(∫_W)((x^2)+((y )^2))dV, where W is the solid bounded between the two paraboloids z=(x^2)+(y^2) and z=8-(x^2)-(y^2).

8. Evaluate the triple integral (∫_(-2)^2)(∫_(-(Sqrt[4-(x^2)]))^(Sqrt[4-(x^2)]))(∫_(Sqrt[(x^2)+(y^2)])^(Sqrt[4-(x^2)-(y^2)]))((x^2)+(y^2)+(z^2))dzdydx by converting to spherical coordinates.

9. The solid defined by the spherical equation ρ=sin ϕ is called the torus.
a. Plot the torus.
b. Calculate the volume of the torus.

10. Ice-Cream Cone: A solid W in the shape of an ice-cream cone is bounded below by the cylinder z=(Sqrt[(x^2)+(y^2)]) and above by the sphere (x^2)+(y^2)+(z^2)=8.  Plot W and determine its volume.

15.5  Applications of Multiple Integrals

Students should read Section 15.5 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Mass as Double Integral: Consider a lamina (thin plate) D in  (ℝ^2) with continous mass density ρ(x,y).  Then the mass of D is given by the double integral

M=∫(∫_D)ρ(x,y)dA

where the domain of integration is given by the region that describes the lamina D.

Example 15.11.  Calculate the mass of the lamina D bounded between the parabola y=(x^2) and y=4 with density ρ(x,y)=y.

Solution: Here is a plot of the lamina D (shaded):

Plot[{x^2,4},{x,-2,2},ImageSize->{250},Filling->{2->{1}}]

We can view D as a Type I region described by -2≤x≤2 and (x^2)≤y≤4.  Thus. the mass of the lamina is given by the double integral:

Integrate[y,{x,-2,2},{y,x^2,4}]

NOTE: Mass of a lamina can also be interpreted as the volune of the solid bounded by its density function over D as shown in the following plot:

Plot3D[y,{x,-2,2},{y,0,4},RegionFunction->(#1^2<#2<4&),Filling->Bottom,Mesh->None,ImageSize->250]

Example 15.12.  Suppose a circular metal plate D, bounded by  (x^2)+(y^2)=9, has electrical charge density ρ(x,y)=(Sqrt[9-(x^2)-(y^2)]).  Calculate the total charge of the plate.

Solution: Here is a plot of the metal plate D (shaded):

Integrate[y,{x,-2,2},{y,x^2,4}]

Plot[{Sqrt[9-x^2],-Sqrt[9-x^2]},{x,-3,3},ImageSize->{250},Filling->{2->{1}},AspectRatio->1]

We shall calculate the total charge of the plate using polar coordinates, which will simplify the corresponding double integral.  Since ρ(r,θ)=(Sqrt[9-(r^2)])and D is a simple polar region described by r=3, the total charge is

Integrate[Sqrt[9-r^2]*r,{r,0,3},{theta,0,2Pi}]

Mass as Triple Integral: We can extend the notion of mass to a solid region W in (ℝ^3).  Suppose W is bounded between two surfaces z=ψ(x,y) and z=ϕ(x,y), where both are defined on the same domain D with ψ(x,y)≤ϕ(x,y), and has density ρ(x,y,z).  Then the mass of W can be expressed by the triple integral

M=∫∫(∫_W)ρ(x,y,z)dV=∫(∫_D)(∫_(ψ(x,y))^(ϕ(x,y)))ρ(x,y,z)dzdA

Example 15.13.  Calculate the mass of the solid region W bounded between the planes z=1-x-y and z=1+x+y and situated over the triangular domain D bounded by x=0, y=0, and y=1-x.  Assume the density of W is given by ρ(x,y,z)=1+(x^2)+(y^2).

Solution: Here is a plot of the solid region W:

Plot3D[{1-x-y,1+x+y},{x,0,1},{y,0,1-x},ViewPoint->{1,1,1},Filling->{1->1,2->1},Ticks->{Automatic,Automatic,{1,2}},ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

The mass of the solid is given by the triple iterated integral (∫_0^1)(∫_0^(1-x))(∫_(1-x-y)^(1+x+y))(1+(x^2)+(y^2))dzdydx:

Integrate[1+x^2+y^2,{x,0,1},{y,0,1-x},{z,1-x-y,1+x+y}]

Center of Mass: Given a lamina D in (ℝ^2), its center of mass ((Subscript[x, CM]),(Subscript[y, CM])) (or balance point) is defined as the ratio of its moments (with respect to the coordinate axes) to its mass:

(Subscript[x, CM])=((Subscript[M, y])/M), (Subscript[y, CM])=((Subscript[M, x])/M)

where the moments (Subscript[M, y]) and (Subscript[M, x]) are defined by

(Subscript[M, y])=(1/A)∫(∫_D)x ρ(x,y)dA,    (Subscript[M, x])=(1/A)∫(∫_D)y ρ(x,y)dA

NOTE: In case the lamina has uniform density, that is, ρ(x,y)=1, then the center of mass is the same as the centroid whose coordinates represent averages of the coordinates over the lamina.

Center of mass (and centroid) can be naturally extended to solid objects in (ℝ^3).  Refer to your textbook for further details.

Example 15.14.  Calculate the mass of the solid region W bounded between the planes z=1-x-y and z=1+x+y and situated over the triangular domain D bounded by x=0, y=0, and y=1-x.  Assume the density of W is given by ρ(x,y,z)=1+(x^2)+(y^2).

Exercises

In Exercises 1 and 2, find the mass of the given lamina D.
1. D is bounded between y=sin (π x) and y=0 along the interval [0,1] and has density ρ(x,y)=x(1-x) .
2. D is bounded by the lines y=x+1, y=-2x-2, and x=1 and has density ρ(x,y)=1+(y^2) .

3. Find the center of mass of the lamina D in Exercises 1 and 2.

4. Find the centroid of the lamina in Exercises 1 and 2.  Compare the centroid of each lamina with its center of mass.

In Exercises 5 and 6, find the mass of the given solidi object W.
5. W is the interior of the tetrahedron enclosed by the planes x=0, y=0, z=0, and z=1-x-y and has density ρ(x,y,z)=1-z.
6. W is the ice-cream cone bounded below by the cylinder z=(Sqrt[(x^2)+(y^2)]) and above by the sphere (x^2)+(y^2)+(z^2)=8 and has density ρ(x,y,z)=(z^2).

7. Find the center of mass of the tetrahedron in Exercises 5 and 6.  Refer to your textbook for appropriate formulas.

8. Find the centroid of the tetrahedron in Exercises 5 and 6.  Compare this with its center of mass.  Refer to your textbook for appropriate formulas.

15.6  Change of Variables

Students should read Section 15.6 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

A change of variables is often useful for simplifying integrals of a single variable (commonly referred to as u-substitution):

(∫_a^b)f(x)dx=(∫_c^d)f(g(u))g'(u)du

where x=g(u), a=g(c), and b=g(d).  This substitution formula allows one to transformation an integral in the variable x to one in a new variable u.  Observe that the interval [c,d] is mapped to interval [a,b] under the function g.

This technique can be extended to double integrals of the form ∫(∫_D)f(x,y)dxdy, where a change of variables is described by a transformation G(u,v)=(x,y), which maps a region (Subscript[D, 0]) in the uv-coordinate plane to the region D in the xy-coordinate plane.

The following Change of Variables Formula converts a double integral from the xy-coordinate system to a new coordinate system defined by u and v:

Change of Variables Formula (Coordinate Transformation): If G(u,v)=(x(u,v),y(u,v)) is a (C^1)-mapping from (Subscript[D, 0]) to D, then

∫(∫_D)f(x,y)dxdy=∫(∫_(Subscript[D, 0]))f(x(u,v),y(u,v))|((∂(x,y))/(∂(u,v)))|dudv

where |((∂(x,y))/(∂(u,v)))|, referred to as the Jacobian of G and also denoted by Jac(G), is given by

Jac(G)=|((∂(x,y))/(∂(u,v)))|=|({{((∂x)/(∂u)), ((∂x)/(∂v))}, {((∂y)/(∂u)), ((∂y)/(∂v))}})|=((∂x)/(∂u))((∂y)/(∂v))-((∂x)/(∂v))((∂y)/(∂u))

The Jacobian relates the area of any infinitesimal region inside (Subscript[D, 0]) with the corresponding region inside D=G((Subscript[D, 0])).  In fact, if G is a linear map, then Jac(G) is constant and is equal in magnitude to the ratio of the areas of D to that of (Subscript[D, 0]):

Jacobian of a Linear Map: If G(u,v)=(A u+C v,B u+D v) is a linear mapping from (Subscript[D, 0]) to D, then Jac(G) is constant with value

Jac(G)=|({{A, C}, {B, D}})|=A D-B C

Moreover,

Area(D)=|Jac(G)|Area((Subscript[D, 0]))

Refer to your textbook for a detailed discussion of transformations of plane regions.

Example 15.12.  Make an appropriate changes of variables to calculate the double integral ∫(∫_D)x ydA, where D is the region bounded by the curves x y=1, x y=2, x (y^2)=1, and x (y^2)=2.

Solution: Here is a plot of the shaded region D bounded by the four given curves:

plot1=ContourPlot[{x*y==1,x*y==2,x*y^2==1,x*y^2==2},{x,0,5},{y,0,5},AspectRatio->Automatic,ImageSize->{250}];
plot2=ContourPlot[1,{x,0,5},{y,0,5},AspectRatio->Automatic,RegionFunction->Function[{x,y},1<x y<2 &&1<x y^2<2],ImageSize->{250},PlotPoints->100];
Show[plot1,plot2]

Observe that D is rather complicated.  Since D can be described by the inequalities 1<x y<2 and 1<x (y^2)<2, we make the natural change of variables u=x y and v=x (y^2), which transforms D to a simple square region (Subscript[D, 0]) in the uv-plane bounded by u=1, u=2, v=1, and v=2:

ContourPlot[1,{u,0,3},{v,0,3},ImageSize->{250},RegionFunction->Function[{u,v},1<u<2 &&1<v<2]]

To find the formula for our transformation G(u,v)=(x(u,v),y(u,v)) that maps (Subscript[D, 0]) to D, we solve for x and y in terms of u and v:

Clear[sol,x,y,u,v]
sol=Solve[{u==x*y,v==x*y^2},{x,y}]

It follows that G(u,v)=((u^2)/v,v/u) and the corresponding Jacobian is

x=sol[[1,1,2]]
y=sol[[1,2,2]]
Jac=D[x,u]*D[y,v]-D[x,v]*D[y,u]

Thus, the given integral transforms to ∫(∫_D)x ydA=∫(∫_(Subscript[D, 0]))(u/v)dA=(∫_1^2)(∫_1^2)(u/v)dvdu with value

Integrate[u/v,{u,1,2},{v,1,2}]

Exercises

1. Consider the transformation G(u,v)=(2u+v,u-3v).
a. Set D=G((Subscript[D, 0])) where (Subscript[D, 0])={0≤u≤1,0≤v≤2}.  Make a plot of D and describe its shape.
b. Compute Jac(G).
c. Compare the area of D with that of (Subscript[D, 0]).  How does this relate to Jac(G)?

2. Compute the area of the ellipse ((x^2)/4)+((y^2)/9)=1 by viewing it as a transformation of the unit circle (u^2)+(v^2)=1 under a linear map G(u,v)=(x(u,v),y(u,v)) and using the area relationship described by Jac(G).

3. Evaluate the integral ∫(∫_D)x ydA, where D is the region in the first quadrant bounded by the equations y=x, y=4x, x y=1, and x y=4.  HINT: Consider the change of variables u=x y and v=y.

4. Evaluate the integral ∫(∫_D)(x+y)/(x-y)dA, where D is the parallelogram bounded by the lines x-y=1, x-y=3, 2x+y=0, and 2x+y=2.  HINT: Consider the change of variables u=x-y and v=2x+y.

5. Evaluate the integral ∫(∫_D)(y/x)dA, where D is the region bounded by the circles (x^2)+(y^2)=1, (x^2)+(y^2)=4 and lines y=x, y=3x.  HINT: Consider the change of variables u=(x^2)+(y^2) and v=y/x.

Chapter 16	Line and Surface Integrals

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

16.1  Vector Fields

Students should read Section 16.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let (Subscript[F, 1]), (Subscript[F, 2]), and (Subscript[F, 3]) be functions of x, y, and z. The vector-valued function

F(x,y,z)=⟨(Subscript[F, 1])(x,y,z),(Subscript[F, 2])(x,y,z),(Subscript[F, 3])(x,y,z)⟩

is called a vector field. We have already encountered a vector field in the form of the gradient of a function. Other useful examples of vector fields are the gravitational force, the velocity of fluid, magnetic fields, and electric fields.

We use the Mathematica commands VectorFieldPlot and VectorFieldPlot3D to plot the graphs of vector fields.  However, before using these commands, it is advisable to load the VectorFieldPlots package. This is done by evaluating

Needs["VectorFieldPlots`"]

Example 16.1.  Draw the following vector fields.
a) F(x,y)=⟨sin y,cos x⟩ 	b)   F(x,y,z)=⟨y, x+z,2x-y⟩

Solution:
a)

Clear[F,x,y,z]
 F[x_,y_]={Sin[y],Cos[x]}

VectorFieldPlot[F[x,y],{x, -5,5},{y,-4,4},ImageSize->{250} ]

Here is another display of the preceding vector field with some options specified.

VectorFieldPlot[F[x,y],{x, -5,5},{y,-4,4},Axes->True,AxesOrigin ->{0,0},Frame->False,ColorFunction->Hue,ImageSize->{250}]

To see other available options of VectorFieldPlot, evaluate the command Options[VectorFieldPlot].

b) We shall use two of the options of VectorFieldPlot3D, which does not have as many options as VectorFieldPlot. (Again, you can find these by evaluating Options[VectorFieldPlot3D].)

Clear[F,x,y,z]
F[x_,y_,z_]={y (z^2),x (z^2),2 x y z}
VectorFieldPlot3D[F[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},ColorFunction->Hue,VectorHeads->True,ImageSize->{250}]

Example 16.2.  Draw the unit radial vector fields:
a) F(x,y)=⟨(x/(Sqrt[(x^2)+(y^2)])),(y/(Sqrt[(x^2)+(y^2)]))⟩	b)   F(x,y,z)=⟨(x/(Sqrt[(x^2)+(y^2)+(z^2)])),(y/(Sqrt[(x^2)+(y^2)+(z^2)])),(z/(Sqrt[(x^2)+(y^2)+(z^2)]))⟩

Solution: For convenience, we define both vector fields to be 0 at the origin. We shall use the If command to do so.

a)

Clear[F,x,y]
F[x_,y_]=If[(x^2)+(y^2)≠0,(({x,y})/(Sqrt[(x^2)+(y^2)])),{0,0}]
VectorFieldPlot[F[x,y],{x,-3,3},{y,-3,3},ImageSize->{250}]

b)

Clear[F,x,y,z]
F[x_,y_,z_]=If[(x^2)+(y^2)+(z^2)≠0,(({x,y,z})/(Sqrt[(x^2)+(y^2)+(z^2)])),{0,0,0}]
VectorFieldPlot3D[F[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3},ColorFunction->Hue,VectorHeads->True,ImageSize->{250}]

Exercises

In Exercise 1 through 4, draw the given vector fields.
1.   F(x,y)=⟨ (y^2)-2xy, xy+6(x^2)⟩			2.   F(x,y,z)=⟨sin x,cos y ,xz⟩   
3.   F(x,y)=⟨  -(y/(Sqrt[(x^2)+(y^2)])),(x/(Sqrt[(x^2)+(y^2)]))⟩			4,   F(x,y,z)=⟨ x+cos (xz),y sin (xy),xz cos (yz) ⟩

In Exerices 5 and 6, calculate and plot the gradient vector field for each of the following functions.
5.   f(x,y)=ln(x+(y^2))				 6.  f(x,y,z)=sin x  (cos z/y)

16.2  Line Integrals

Students should read Section 16.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Suppose C is a smooth curve in space whose parametric equations are given by

x=x(t), y=y(t),  z=z(t)

where a≤t≤b.  Let (Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3]), .... , (Subscript[C, N]) be a partition of the curve C with arc length (Subscript[Δs, 1]), (Subscript[Δs, 2]), (Subscript[Δs, 3]), ...  , (Subscript[Δs, N]) and let (Subscript[P, 1]), (Subscript[P, 2]), (Subscript[P, 3]), ... , (Subscript[P, N]) be points on the subarcs.

If f(x,y,z) is a function that is continuous on the curve C, then the line integral of f is defined by

(∫_C)f(x,y,z)ds=(lim_((Subscript[Δs, i])->0))(∑_(i=1)^N)f((Subscript[P, i]))(Subscript[Δs, i])

NOTE: If c(t)=⟨x(t), y(t),z(t)⟩ is the vector equation of the curve C, then it can be shown (refer to your calculus textbook) that

(∫_C)f(x,y,z)ds=(∫_a^b)f(c(t)‖c'(t)‖dt

In addition, if F(x,y,z)=⟨(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])⟩ is a vector field that is continuous on C, then the line integral of F over C is given by

(∫_C)F(x,y,z)·ds=(∫_C)(F·T)ds=(∫_a^b)F(c(t))·c'(t)dt

where T is the unit vector T=((c'(t))/(‖c'(t)‖)) and F·T is the dot product of F and T.

Example 16.3.  Find  (∫_C)f(x,y,z)ds, where f(x,y,z)=x y+(z^2) and C is given by  x= t, y= (t^2), and z=(t^3), for 0≤t≤1.

Solution:

Clear[x,y,z,t,f,c]
f[x_,y_,z_]=((x )^2)y+x z
x[t_]=t
y[t_]=(t^2)
z[t_]=(t^3)
c[t_]={x[t],y[t],z[t]}

(∫_0^1)f[x[t],y[t],z[t]]Norm[c'[t]]\[DifferentialD]t

Here is a numerical approximation of the preceding line integral.

NIntegrate[f[x[t],y[t],z[t]]Norm[c'[t]],{t,0,1}]

Example 16.4.  Find  (∫_C)F(x,y,z)·ds, where F(x,y,z)=⟨ x  z, z (y^2), y (x^2)⟩ and the curve C is given by  x= 2t, y= sin t, and z=cos t, 0≤t≤2π.

Solution:

Clear[x,y,z,t,f,c]
F[x_,y_,z_]={x z,z (y^2), y (x^2)} 
x[t_]=2t 
y[t_]= Sin[t]
z[t_]=Cos[t]
c[t_]={x[t],y[t],z[t]}

(∫_0^(2Pi))F[x[t],y[t],z[t]].c'[t]\[DifferentialD]t

N[%]

Exercises

1.  Find  (∫_C)f(x,y,z)ds, where:
a.  f(x,y,z)= x (y^(2 ))-4 zy and C is given by  x= 2t, y= (t^(2/3)), and z=1-3(t^2), for 0≤t≤1.
b.  f(x,y,z)=(yz/x) and C is given by  x= ln t, y= (t^2), and z=3t, for 3≤t≤5.

2.   Find  (∫_C)F(x,y)·ds, where:
a.   F(x,y)=⟨ (e^(3x-2y)),(e^(2x+3y))⟩ and C is given by  x= 2t, y= sin t,  0≤t≤ π
b.  F(x,y)=⟨(x^2),yx+(y^2)⟩ and C is the unit circle center at the origin.

3.  Find  (∫_C)F(x,y,z)·ds, where:
a.   F(x,y,z)=⟨xyz,-xz,xy ⟩ and C is given by  x= t, y= 2(t^2), z=3t  0≤t≤ 1
b.   F(x,y,z)=⟨(xy^3),z+(x^2),(z^3)⟩ and C is the line segment joining (-1,2,-1) and (1, 3, 4).

16.3  Conservative Vector Fields

Students should read Section 16.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let F(x,y,z)=⟨(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])⟩ be a vector field. Let (Subscript[C, 1]) and (Subscript[C, 2]) be any two different curves with the same initial point P and end point Q.  We say that the vector field F is path independent if

(Subscript[(∫_C), 1])F(x,y,z)·ds=(Subscript[(∫_C), 2])F(x,y,z)·ds

A vector field that is path independent is called conservative.

NOTE 1: A vector field F is conservative if

(∫_C)F(x,y,z)·ds=0

for every closed curve C.

NOTE 2: If F=∇u is the gradient of a function u=u(x,y,z), then we say that u is the potential of F.  Moreover, if the end points of C are P and Q, we have

(∫_C)F(x,y,z)·ds=u(P)-u(Q)

In particular, if the curve is closed, that is, if P=Q, then

(∫_C)F(x,y,z)·ds=0

Therefore, gradient is conservative. The converse of this statement is true if its domain is an open connected domain.

NOTE 3: Let F=⟨(Subscript[F, 1]),(Subscript[F, 2])⟩.  If F=∇u=⟨((∂ u)/(∂ x)),((∂ u)/(∂ y))⟩, then (Subscript[F, 1])=((∂ u)/(∂ x)) and (Subscript[F, 2])=((∂ u)/(∂ y)). Taking the partial derivative of (Subscript[F, 1])with respect to y and that of (Subscript[F, 2]) with respect to x and using the fact that (((∂^2) u)/(∂ x∂ y))=(((∂^2) u)/(∂ y∂ x)), we see that (Subscript[F, 1]) and (Subscript[F, 2]) must satisfy

((∂ (Subscript[F, 1]))/(∂ y))=((∂ (Subscript[F, 2]))/(∂ x))

This equation is used to check if a vector field is conservative.  In that case, we solve (Subscript[F, 1])=((∂ u)/(∂ x)) for u by integrating with respect to x and then use the equation (Subscript[F, 2])=((∂ u)/(∂ y)) to find the constant of integration. Here is an example.

Example 16.5.  Show that the vector function F=⟨3(x^2)-2xy +2,6(y^2)-(x^2)+3⟩  is conservative and find its potential.

Solution: Here, (Subscript[F, 1])=x (y^2) and (Subscript[F, 2])=(x^2)y. We now compare ((∂ (Subscript[F, 1]))/(∂ y)) and ((∂ (Subscript[F, 2]))/(∂ x)) to verify if F is conservative.

Clear[x,y,F1,F2]
F1[x_,y_]= 3(x^2)-2x y+2
F2[x_,y_]=6(y^2)-(x^2)+3

D[F1[x,y],y]
D[F2[x,y],x]

Thus, the vector field is conservative. To find its potential u, we integrate (Subscript[F, 1])=((∂ u)/(∂ x)) with respect to x to get

Clear[h,u]
u=Integrate[F1[x,y],x] +h[y]

Note that the addition of h(y) is necessary because the constant of integration may depend on y.  We now solve the equation (Subscript[F, 2])=((∂ u)/(∂ y)) for h'(y).

Clear[sol]
sol=Solve[D[u,y]==F2[x,y],h'[y]]

This means that h'(y)=3(1+2(y^2)).

Integrate[sol[[1,1,2]],y]

Hence, h(y)=3y+2(y^2)  and so u(x,y)=2x+(x^3)-(x^2) y +3y+2(y^3) is the potential of F.

NOTE 4: Let F=⟨(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])⟩.  If F=∇u=⟨((∂ u)/(∂ x)),((∂ u)/(∂ y)),((∂ u)/(∂ z))⟩, then (Subscript[F, 1])=((∂ u)/(∂ x)), (Subscript[F, 2])=((∂ u)/(∂ y)) and (Subscript[F, 3])=((∂ u)/(∂ z)). Taking the partial derivative of (Subscript[F, 1])with respect to y and that of (Subscript[F, 2]) with respect to x and using the fact that (((∂^2) u)/(∂x∂y))=(((∂^2) u)/(∂y∂x)), we see that (Subscript[F, 1]) and (Subscript[F, 2]) must satisfy

((∂ (Subscript[F, 1]))/(∂ y))=((∂ (Subscript[F, 2]))/(∂ x))

Taking the partial derivative of (Subscript[F, 1])with respect to z and that of (Subscript[F, 3]) with respect to x and using the fact that (((∂^2) u)/(∂x∂z))=(((∂^2) u)/(∂z∂x)), we see that (Subscript[F, 1]) and (Subscript[F, 3]) must satisfy

((∂ (Subscript[F, 1]))/(∂ z))=((∂ (Subscript[F, 3]))/(∂ x))

The preceding two equations can be used to check if a vector field is conservative. If this the case, we solve (Subscript[F, 1])=((∂ u)/(∂ x)) for u by integrating with respect to x and then use (Subscript[F, 2])=((∂ u)/(∂ y)) to find the constant of integration. We show this by the following example.

Example 16.6.  Show that the vector function F=⟨yz+yz cos (xy ), xz+xz cos (xy) ,xy+sin (xy)⟩  is conservative and find its potential.

Solution: Here, (Subscript[F, 1])=y z+y z cos( x y) ,  (Subscript[F, 2])=x z +x z cos (x y),  and (Subscript[F, 3])=x y + sin(x y) .

Clear[x,y,F1,F2,F3]
F1[x_,y_,z_]= y z +y z Cos[x y]
F2[x_,y_,z_]= x z + x z Cos[x y]
F3[x_,y_,z_]= x y +  Sin[x y]

We now compare ((∂ (Subscript[F, 1]))/(∂ y)) and ((∂ (Subscript[F, 2]))/(∂ x)):

D[F1[x,y,z],y]
D[F2[x,y,z],x]

Next, we compare ((∂ (Subscript[F, 1]))/(∂ z)) and ((∂ (Subscript[F, 2]))/(∂ x)):

D[F1[x,y,z],z]
D[F3[x,y,z],x]

Thus, the vector field is conservative. To find its potential u, we integrate (Subscript[F, 1])=((∂ u)/(∂ x)) with repsetct to x to get

Clear[u, h]
u=Integrate[F1[x,y,z],x] +h[y,z]

Note that the addition of h(y,z) is necessary because the constant of intgeration can depend on y and z. We now solve the equation (Subscript[F, 2])=((∂ u)/(∂ y)) for ((∂ h)/(∂ y)).

Clear[sol]
sol=Solve[D[u,y]==F2[x,y,z],(∂_y) h[y,z]]

This means that ((∂ h)/(∂ y))=0 and hence h is a function of z only.   Next, we solve the equation (Subscript[F, 3])=((∂ u)/(∂ z)) for ((∂ h)/(∂ z)).

Clear[sol2]
sol2=Solve[D[u,z]==F3[x,y,z],(∂_z) h[y,z]]

Hence, ((∂ h)/(∂ z))=0 and we can take h=0. Therefore, u=x y z +z sin (x y) is the potential for the vector field F.

Exercises

1. Show that the vector field F=⟨ (y^3)-3(x^2)y, 3(xy^2)-(x^3)⟩  is conservative and find its potential.

2. Show that the vector field F=⟨y z+((2 xy)/z), xz+((x^2)/z), xy-(((x^2)y)/(z^2))⟩  is conservative and find its potential.

3. Determine whether the vector field F=⟨(x^2), yx +(e^z), y (e^z)⟩  is conservative.  If it is, find its potential.

16.4  Parametrized Surfaces and Surface Integrals

Students should read Section 16.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

A parametrized surface is a surface whose points are given in the form

G(u,v)=(x(u,v), y(u,v), z(u,v))

where u and v (called parameters) are independent variables used to describe a domain D (called the parameter domain).

The command for plotting parametrized surfaces is ParametricPlot3D.  This command has been discussed in Section 14.1.2 of this text.

Example 16.7.  Plot the parametrized surface defined by G(u,v)=(cos u sin v, 4 sin u cos v, cos v) over the domain D={(u,v)| 0≤u≤2π, 0≤v≤2π}.

Solution:

ParametricPlot3D[{ Cos[u] Sin[v],4 Sin[u] Cos[v], Cos[v]},{u,0,2Pi},{v,0,2Pi},ImageSize->{250}]

Example 16.8.  Plot the parametrized surface defined by G(u,v)=(u cos v,u sin v,1-(u^2)) over the domain D={(u,v)|0≤u≤1,0≤v≤2π}.

Solution:

ParametricPlot3D[{u Cos[v]  ,u Sin[v]  , 1-(u^2)},{u,0,1},{v,0,2Pi},ColorFunction->"BlueGreenYellow",ImageSize->{250},
ImagePadding->{{15,15},{15,15}}]

NOTE: On a parametrized surface G(u,v)=(x(u,v), y(u,v), z(u,v)), if we fix one of the variables, we get a curve on the surface.  The plot following shows the curves corresponding to u=3/4 (latitude) and v=5π/3 (longitude).

Clear[plot1,plot2,plot3] 
plot1= ParametricPlot3D[{u  Cos[v]  ,u Sin[v]  , 1-(u^2)},{u,0,1},{v,0,2Pi},ColorFunction->"BlueGreenYellow"];
plot2=ParametricPlot3D[{  3/4 Cos[v]  ,3/4  Sin[v]  , 7/16},{v,0,2Pi},PlotStyle->{Thickness[0.01],Red}];
plot3=ParametricPlot3D[{  u  Cos[5Pi/3]  ,u  Sin[5Pi/3]  ,1-(u^2)},{u,0,1},PlotStyle->{Thickness[0.01],Blue}];
Show[plot1,plot2,plot3,PlotRange->All,ImageSize->{250},
ImagePadding->{{15,15},{15,15}}]

Let P=G((Subscript[u, 0]),(Subscript[v, 0])) be a point on the parametrized surface S.  For fixed v=(Subscript[v, 0]), the tangent vector to the curve G(u,(Subscript[v, 0])) at ((Subscript[u, 0]),(Subscript[v, 0]))  is given by

(Subscript[T, u])=((∂G)/(∂u))((Subscript[u, 0]),(Subscript[v, 0]))

while the tangent vector for  G((Subscript[u, 0]),v) corresponding to a fixed u=(Subscript[u, 0]) is given by

(Subscript[T, v])=((∂G)/(∂v))((Subscript[u, 0]),(Subscript[v, 0]))

These two vectors are tangent to the surface S.  Thus, the normal vector n to the tangent plane at G((Subscript[u, 0]),(Subscript[v, 0])) is given by

n(P)=n((Subscript[u, 0]),(Subscript[v, 0]))=(Subscript[T, u])×(Subscript[T, v])

Example 16.9.  Consider the parametrized surface G(u,v)=( u cos v, u sin v, 1-(v^2)).
a) Find (Subscript[T, u]), (Subscript[T, v]), and n.
b) Find the equation of the tangent plane at (1/2, 5π/3).
c) Plot the tangent plane and surface.

Solution:  Let us define G as a function of u and v in Mathematica.

Clear[G,u,v]
G[u_,v_]={u Cos[v], u Sin[v], 1-(u^2)}

a) We use Tu for (Subscript[T, u]) and Tv for (Subscript[T, v]). We evaluate these as functions of u and v.

Clear[Tu,Tv,n]
Tu[u_,v_]=D[G[u,v],u]
Tv[u_,v_]=D[G[u,v],v]
n[u_,v_]=Cross[Tu[u,v],Tv[u,v]]

b) The normal vector to the tangent plane at (1/2, 5π/3) is

Clear[normal]
normal=n[1/2,5Pi/3]

The tangent plane passes through the point

Clear[point]
point=G[1/2,5Pi/3]

Thus, the equation of the tangent plane is given by

Clear[tplane]
tplane=normal.({x,y,z}-point)==0

which simplifies to

Simplify[tplane]

c) Here is the plot of the surface and the tangent plane. Observe that we have used ColorFunction and ColorFunctionScaling options.

Clear[plot1, plot2]
plot1= ParametricPlot3D[G[u,v],{u,0,1},{v,0,2Pi},ColorFunction->"BlueGreenYellow"];
plot2=ContourPlot3D[2 x+4 z==5+2 (Sqrt[3]) y, {x, -3,3}, {y, -3,3},{z, -4,4},ColorFunction->Function[{x,y,z},Hue[Mod[z,1]]],ColorFunctionScaling->False];
Show[plot1, plot2,ImageSize->{250},ImagePadding->{{15,15},{15,15}}]

NOTE: The area A(S) of a parametrized surface S: G(u,v)=(x(u,v), y(u,v), z(u,v)), where (u,v)∈D, is given by

A(S)=∫(∫_D)||n(u,v)||dudv

If f(x,y,z) is continuous at all points of S, then the surface area of f over S is given by

∫(∫_S)f(x,y,z)dS=∫(∫_D)f(G(u,v))||n(u,v)||dudv

Example 16.10.  Show the following:
a) The area of the cylinder of height h  and radius r is 2π rh.
b) The area of the sphere of radius r is 4π (r^2).

Solution: 
a) A parametric equation of the cylinder of height h  and radius r can be given by

x=r cos v, y=r sin v, and z=u, where  0≤v≤2π , 0≤u≤h

Thus, the cylinder is given by G(u,v)=(r cos u, r sin u, v).

Clear[ G,u,v,r]
G[u_,v_]={r Cos[v], r Sin[v],u}

Here is a plot of the cylinder with r=3 and h=5:

r=3;h=5;
ParametricPlot3D[G[u,v], {u, 0, h},{v, 0, 2Pi}]

To compute the surface area of the cylinder, we need to compute its normal vector.

Clear[Tu,Tv,n,r,h]
Tu[u_,v_]=D[G[u,v],u];
Tv[u_,v_]=D[G[u,v],v];
n[u_,v_]=Cross[Tu[u,v],Tv[u,v]]

Here is a plot of the cylinder with its normal vector for  r=3 and h=5:

r=3;h=5;
Clear[plot1,plot2]
plot1=ParametricPlot3D[G[u,v], {u, 0, h},{v, 0, 2Pi}];
plot2=VectorFieldPlot3D[n[u,v],{u,0,h},{v,-2Pi,2Pi},{z,-3,3},  VectorHeads->True,PlotPoints->15];
Show[plot1,plot2,ImageSize->{250}]
Clear[r,h]

The surface area is

SArea=(∫_0^h)(∫_0^(2Pi))Norm[n[u,v]]\[DifferentialD]v\[DifferentialD]u

Since r>0, |r|=r and hence the preceding output is 2π r h .

b) A parametric equation of the sphere of radius r is

x=r cos u sin v, y=r sin u sin v, z=r cos v

where 0≤u≤2π and 0≤v≤π.  Thus, the sphere is given by G(u,v)=(r cos u sin v,r sin u sin v,r cos v).

Clear[ G,u,v,r]
G[u_,v_]={r Cos[u] Sin[v], r Sin[u] Sin[v],r Cos[v]}

Here is a plot of the sphere with r=3.

r=3; 
ParametricPlot3D[G[u,v],{u, 0, 2Pi}, {v, 0, Pi},ImageSize->{250}]

To compute the surface area of the sphere, we need to compute its normal vector.

Clear[Tu,Tv,n,r]
Tu[u_,v_]=D[G[u,v],u];
Tv[u_,v_]=D[G[u,v],v];
n[u_,v_]=Cross[Tu[u,v],Tv[u,v]]

Here is a plot of the sphere with its normal vector for  r=3.

r=3;h=5;
Clear[plot1,plot2]
plot1=ParametricPlot3D[G[u,v],{u, 0, 2Pi}, {v, 0, h}];
plot2=VectorFieldPlot3D[n[u,v],{u,-2Pi,2Pi},{v,0,h},{z,-3,3},  VectorHeads->True,PlotPoints->10];
Show[plot1,plot2,ImageSize->{250}]
Clear[r,h]

The surface area is

SArea=(∫_0^Pi)(∫_0^(2Pi))Norm[n[u,v]]\[DifferentialD]u\[DifferentialD]v

For a real number r, the conjugate of r is r and hence the preceding output is 4π (r^2).

Example 16.11.  Consider the parametrized surface S defined by G(u,v)=( u cos v, u sin v, v), where 0≤u≤1, 0≤v≤2π.
a) Find the surface area of S.
b) Evaluate ∫(∫_S)xyzdS.

Solution:
a)

Clear[ G,u,v]
G[u_,v_]={u Cos[v], u Sin[v],v}

Clear[Tu,Tv,n]
Tu[u_,v_]=D[G[u,v],u]
Tv[u_,v_]=D[G[u,v],v]
n[u_,v_]=Cross[Tu[u,v],Tv[u,v]]

The surface area A(S) is given by

SArea=(∫_0^1)(∫_0^(2Pi))Norm[n[u,v]]\[DifferentialD]v\[DifferentialD]u

which is approximately equal to

N[%]

b) We define f:

Clear[f]
f[x_,y_,z_]=x y z

The surface integral of f is

(∫_0^1)(∫_0^(2Pi))f[G[u,v][[1]],G[u,v][[2]],G[u,v][[3]]]Norm[n[u,v]]\[DifferentialD]v\[DifferentialD]u

Or numerically,

N[%]

Exercises

1. Plot the parametrized surface G(u,v)=((e^u)sin v,(e^u) cos v,v) over the domain D={(u,v)| -1≤u≤1, 0≤v≤2π}.

2. Plot the parametrized surface G(u,v)=(3sin u cos v,sin u sin v,cos v+3 cos u) over the domain D={(u,v)| 0≤u≤2π, 0≤v≤2π}.

3. Consider the parametrized surface G(u,v)=((e^(-u))  cos v,(e^u)sin v,(e^u)cos v).
a.  Find (Subscript[T, u]), (Subscript[T, v]), and n.
b.  Find the equation of the tangent plane at (0, π/ 2).
c.  Plot the tangent plane and surface.

4. Consider the parametrized surface S: G(u,v)=( u-v ,3u+v,(u^2)-2 u v +6(v^2)), where 0≤u≤1, 0≤v≤1.
a.  Find the surface area of S. (Use NIntegrate for faster integration.)
b.  Evaluate ∫(∫_S)(3x+ 2(y^2)-  (z^2))dS.

16.5  Surface Integrals of Vector Fields

Students should read Section 16.5 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

An orientation of a surface S is a continuously varying choice of the unit normal vector (Subscript[e, n])(P) at each point of the surface.  Thus, (Subscript[e, n]) is given by either

(Subscript[e, n])(P)=((n(P))/(||n(P)||))   or   (Subscript[e, n])(P)=-((n(P))/(||n(P)||))

If F(x,y,z) is continuous at all points of a parametrized surface S, then the surface integral of F over S is given by

∫(∫_S)F·dS=∫(∫_S)(F·(Subscript[e, n]))dS

where (Subscript[e, n]) is the unit normal determined by an orientation. The surface integral of F is also called the flux of F across S.

The surface integral of F over a parametrized surface S given by G(u,v)=(x(u,v),y(u,v),z(u,v)), where (u,v)∈D, is given by

∫(∫_S)F·dS=∫(∫_S)(F·(Subscript[e, n]))dS =∫(∫_D)F(G(u,v))·n(u,v) dudv

Example 16.12.  Find  ∫(∫_S)F·dS, where F(x,y,z)=⟨x  z, z , y x⟩ and S is given by  G(u,v)=(u-(v^2),u v,(u^2)-v), 0≤u≤2, and 1≤v≤3.

Solution:

Clear[F,G,x,y,z,u,v]
F[x_,y_,z_]={x z, z, y x}
G[u_,v_]={u-(v^2), u v, (u^2)-v}

Clear[Tu,Tv,n]
Tu[u_,v_]=D[G[u,v],u]
Tv[u_,v_]=D[G[u,v],v]
n[u_,v_]=Cross[Tu[u,v],Tv[u,v]]

Flux=(∫_0^2)(∫_1^3)F[G[u,v][[1]],G[u,v][[2]],G[u,v][[3]]] .n[u,v]\[DifferentialD]v\[DifferentialD]u

Example 16.13.  Find  ∫(∫_S)F·dS, where F(x,y,z)=⟨(x^2), (z^2) , y+ (x^2)⟩ and S is the upper hemisphere (x^2)+(y^(2 ))+(z^2)=4 with outward normal orientation.

Solution:  First, we find the parametric equation of the cylinder. This can be given by x=2 cos u sin v, y=2 sin u sin v, and z=2cos v, where 0≤u≤2π and 0≤v≤π/2.

For the hemisphere to have the outward orientation, we note that n=(Subscript[T, v])×(Subscript[T, u]). With this in mind we compute the flux of F across S through the following steps.

Clear[F,G,x,y,z,u,v]
F[x_,y_,z_]={   (x^2), (z^2), (x^2)+y+(z^3)}
G[u_,v_]= { 2 Cos[u]Sin[v]  , 2 Sin[u] Sin[v]   , Cos[v] }

Clear[Tu,Tv,n]
Tu[u_,v_]=D[G[u,v],u]
Tv[u_,v_]=D[G[u,v],v]
n[u_,v_]=Cross[Tv[u,v],Tu[u,v]]

Flux=(∫_0^(Pi/2))(∫_0^(2Pi))F[G[u,v][[1]],G[u,v][[2]],G[u,v][[3]]] .n[u,v]\[DifferentialD]u\[DifferentialD]v

Exercises

1. Find ∫(∫_S)F·dS, where F(x,y,z)=⟨ (e^z),z,y x⟩ and S is given by G(u,v)=(u v,u-v,u),  0≤u≤2, and -1≤v≤1, and oriented by n=(Subscript[T, u])×(Subscript[T, v]).

2. Find  ∫(∫_S)F·dS, where F(x,y,z)= ⟨ z,x ,y ⟩ and S is the portion of the ellipsoid  ((x^2)/16)+((y^(2 ))/9)+((z^2)/4)=1 for which x≤0,   y≤0, and z≤0 with outward normal orientation.

3. Let S be given by G(u,v)=((1+v cos(u/2))cos u,(1+v cos(u/2))sin u , v sin (u/2) ),  0≤u≤2π, and ((-1)/2)≤v≤(1/2).
a. Plot the surface S. (S is an example of a Mobius strip.)
b. Find the surface area of S.
c. Evaluate ∫(∫_S)((x^2)+ 2(y^2)+3  (z^2))dS.
d. Find the intersection points of S and the xy-plane.
e. For each of the points on the intersection of S and the xy-plane, find the normal vector n.
f. Show that n varies continuously but that n(2π,0)=-n(u,0). (This shows that S is not orientable and hence it is impossible to integrate a vector field over S.)

Chapter 17	Fundamental Theorems of Vector Analysis

Useful Tip: If you are reading the electronic version of this publication formatted as a Mathematica Notebook, then it is possible to view 3-D plots generated by Mathematica from different perspectives.  First, place your screen cursor over the plot.  Then drag the mouse while pressing down on the left mouse button to rotate the plot.

NOTE: In order to perform the operations of curl and divergence on vector fields discussed in this section using Mathematica, it is necessary to first load the VectorAnalysis package:

Needs["VectorAnalysis`"]

The Fundamental Theorem of Calculus for functions of a single variable states that the integral of a function f(x) over an interval [a,b] (domain) can be calculated as the difference of its anti-derivative F(x) at the endpoints (boundary) of the interval:

(∫_a^b)f(x)dx=F(b)-F(a)

This integral relationship between domain and boundary can be generalized to vector fields involving the operations of curl and divergence and is made precise by three theorems that will be discussed in this chapter: Green's Theorem, Stoke's Theorem, and Divergence Theorem.

17.1  Green's Theorem

Students should read Section 17.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let  F(x,y)=⟨P(x,y),Q(x,y)⟩ a vector field continuous on an oriented curve C.  Recall that the line integral of F along C is denoted by

(∫_C)F(x,y,z)·ds=(∫_C) Pdx+Qdy

If c(t)=⟨x(t), y(t),z(t)⟩ is the vector equation of the curve C, then

(Subscript[∲, C]) Pdx+Qdy=  (∫_a^b)(P(x(t),y(t)) ((d x)/(d t))+Q(x(t),y(t)) ((d y)/(d t)))d t

The following is a generalization of the Fundamental Theorem of Calculus to two dimensions, which relates a double integral over a region with a corresponding line integral along its boundary.

Green's Theorem: If C is a simple closed curve oriented counterclockwise and D is the region enclosed, and if P and Q are differentiable and have continuous first partial derivatives, then

(Subscript[∲, C]) Pdx+Qdy=∫(∫_D)(((∂ Q)/(∂ x))-((∂ P)/(∂ y))) dA

Refer to your textbook for a detailed discussion and proof of Green's Theorem.

Example 17.1. Compute the line integral (Subscript[∲, C]) (e^(2x+y))dx+(e^(-y))dy, where C is the boundary of the square with vertices (0,0), (1,0), (1,1), (1,0) oriented counterclockwise.

Solution: We will use Green's Theorem.  Thus, we need to verify that the hypotheses of Green's Theorem hold. To this end, we define the function P and Q and compute their partial derivatives.

Clear[x,y,P,Q]
P[x_,y_]= (E^(2x+y))
Q[x_,y_]=(E^(-y))

D[P[x,y],x]
D[P[x,y],y]
D[Q[x,y],x]
D[Q[x,y],y]

The partial derivatives are continuous inside the square and the curve is oriented counterclockwise. Thus, the hypotheses of Green's Theorem are satisfied.  Note that the region D enclosed by C is given by 0≤x≤1 and 0≤y≤1.

(∫_0^1)(∫_0^1)(D[Q[x,y],x]-D[P[x,y],y])\[DifferentialD]y\[DifferentialD]x

N[%]

NOTE: If we were to solve this using the definition of line integral as discussed in Chapter 16 of this text, we would then need to consider four pieces of parametrization of C and then sum the four integrals.   Toward this end, let us use (Subscript[C, 1]) for the lower edge, (Subscript[C, 2]) for the right edge, (Subscript[C, 3]) for the top edge, and (Subscript[C, 4]) for the left edge of the square.  Here are the parametrizations followed by their line integrals.

Clear[x1,x2,x3,x4,y1,y2,y3,y4,t,F,c1,c2,c3,c4]
F[x_,y_]={P[x,y],Q[x,y] } 
x1[t_]=t
y1[t_]= 0
c1[t_]={x1[t],y1[t]}

x2[t_]=1
y2[t_]= t
c2[t_]={x2[t],y2[t]}

x3[t_]=1-t
y3[t_]= 1
c3[t_]={x3[t],y3[t]}

x4[t_]=0
y4[t_]=1- t
c4[t_]={x4[t],y4[t]}

(∫_0^1)F[x1[t],y1[t]].c1'[t]\[DifferentialD]t+(∫_0^1)F[x2[t],y2[t]].c2'[t]\[DifferentialD]t+(∫_0^1)F[x3[t],y3[t]].c3'[t]\[DifferentialD]t+(∫_0^1)F[x4[t],y4[t]].c4'[t]\[DifferentialD]t

N[%]

Exercises

In Exercises 1 through 4, use Green's Theorem to evaluate the given line integral.

1. (Subscript[∲, C]) (y^2)sin x dx+x ydy, where C is the boundary of the triangle with vertices (0,0), (1,0), (1,1), oriented counterclockwise.

2. (Subscript[∲, C]) 2(x^2) y dx+(x^3) dy, where C is the circle (x^2)+(y^2)=4, oriented counterclockwise.

3. (Subscript[∲, C]) ((x^2)+(y^2)) dx+ y (e^x)dy, where C is the boundary of the region bounded between the parabola y=5-(x^2) and the line y=2x-3, oriented clockwise.

4. (Subscript[∲, C])  (x/((x^2)+(y^2)))dx-(y/((x^2)+(y^2))) dy, where C is the boundary of the quarter-annulus situated between the circles (x^2)+(y^2)=1 and (x^2)+(y^2)=9 in the first quadrant (see plot below), oriented counterclockwise.

5. Let F(x,y)=(2x y+(y^3),(x^2)+3x y+2y).  Use Green's Theorem to demonstrate that the line integral (∫_C)F(x,y,z)·ds=0 for every simple closed curve C.  What kind of a vector field do we call F?

17.2  Stokes's Theorem

Students should read Section 17.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let F(x,y,z)=⟨(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])⟩ be a vector field.  The curl of F, denoted by curl(F) or ∇×F, is defined by

curl(F)=∇×F=|({{i, j, k}, {(∂/(∂x)), (∂/(∂y)), (∂/(∂z))}, {(Subscript[F, 1]), (Subscript[F, 2]), (Subscript[F, 3])}})|=⟨((∂(Subscript[F, 3]))/(∂y))-((∂(Subscript[F, 2]))/(∂z)),((∂(Subscript[F, 1]))/(∂z))-((∂(Subscript[F, 3]))/(∂x)),((∂(Subscript[F, 2]))/(∂x))-((∂(Subscript[F, 1]))/(∂y))⟩

Here, we are using the del or symbol ∇ (nabla) to denote the vector operator ∇= ⟨(∂/(∂x)),(∂/(∂y)),(∂/(∂z))⟩.

The Mathematica command for computing the curl of a vector field F is Curl[F,coordsys], where coordsys is the coordinate system of the vector field.  This is demonstrated in the next example.

The following is a generalization of the Fundamental Theorem of Calculus three dimensions, which relates a surface integral involving curl with a corresponding line integral along its boundary.

Stokes's Theorem: If  F(x,y,z) a vector field with continuous partial derivatives and if S is an oriented surface S with boundary ∂S,  then

(∮_(∂S))F·dS=∫(∫_S)curl(F)·dS

If S is closed, then it has no boundary and hence both integrals are equal to 0.

NOTE: Recall that if the surface S is given by G(u,v)=(x(u,v), y(u,v), z(u,v)), where (u,v)∈D, then (∫_S)curl(F)·dS is given by

∫(∫_S)curl(F)·dS=∫(∫_D)curl(F)(G(u,v)) ·n(u,v) dudv

Refer to your textbook for a detailed discussion and proof of Stokes's Theorem.

Example 17.2.  Find the curl of the vector field F(x,y,z)=⟨ x sin (y z),(e^(x /y)) z ,y (x^2)⟩.

Solution: We use the Curl command:

Clear[F,F1,F2,F3,x,y,z]
F1 =x  Sin[y z]
F2  =(E^(x/y))z
F3 =((x )^2) y
F={F1,F2,F3}

Curl[F,Cartesian[x,y,z]]

NOTE: We obtain the same answer for the curl of F using the explicit formula:

curl={(∂_y) F3-(∂_z)F2,(∂_z)F1 -(∂_x)F3,(∂_x)F2-(∂_y)F1}

Or equivalently,

CurlF={D[F3,y]-D[F2,z],D[F3,x]-D[F1,z],D[F2,x]-D[F1,y]}

Example 17.3.  Let f(x,y,z) be a function of three variables with continuous first and second partial derivatives and let F=∇f be the gradient of f.  Find the curl of the vector field F.

Solution:

Clear[f,F1,F2,F3,x,y,z]
F1 =D[f[x,y,z],x]
F2  =D[f[x,y,z],y]
F3 =D[f[x,y,z],z]
F={F1,F2,F3}

Then the curl of F is

Curl[F,Cartesian[x,y,z]]

To see why the curl is zero, let us examine each partial derivative used in computing the curl of F.

D[F3,y]
D[F2,z]

NOTE: Here, (f^((0,1,1)))[x,y,z] stands for the second partial derivative (Subscript[f, yz]). Thus, the two partial derivatives that appear in the x-component of the curl of F are equal and hence their difference is zero. Similarly, we have

D[F3,x]
D[F1,z]

and

D[F2,x]
D[F1,y]

Example 17.4.  Compute  (∮_(∂S))F·dS, where F(x,y,z)=⟨ x y z, z+3x -3y , ((y )^2)x⟩ and S is the upper hemisphere of radius 4.

Solution: Note that ∂S is a circle of radius 4 lying on the xy-plane. Hence, ∂S can be parametrized by the curve c(t)=(x(t),y(t),z(t)) where

x=4 cos t, y=4sin t, z=0, where  0≤t≤2π

We then use this parametrization to evaluate the line integral (∮_(∂S))F·dS=(∫_0^(2π))F(x(t),y(t),z(t))·c'(t)dt:

Clear[F, x,y,z,t,c,curlF]
F[x_,y_,z_]={x y z, z+3x-3y, (y^2) x}
x[t_]=4Cos[t]
y[t_]= 4Sin[t]
z[t_]=0
c[t_]={x[t],y[t],z[t]}

(∫_0^(2Pi)) F[x[t],y[t],z[t]].c'[t]\[DifferentialD]t

Next, we use Stokes's Theorem to obtain the same answer via the corresponding surface integral.  The parametrization of the upper hemisphere of radius 4 is given by S(u,v)={x(u,v),y(u,v),z(u,v)}, where

x=4cos u sin v,  y=4 sin u sin v, and z=4 cos v,    where    0≤u≤2π, 0≤v≤π/2

We now compute the normal of the upper hemisphere:

Clear[S,u,v,Tu,Tv,n] 
S[u_,v_]:= { 4 Cos[u]Sin[v],4Sin[u] Sin[v],4Cos[v] } 
Tu[u_,v_]:=D[S[u,v],u]
Tv[u_,v_]:=D[S[u,v],v]
n[u_,v_]=Cross[Tv[u,v],Tu[u,v]]

The curl of F is

curlF[x_,y_,z_]=Curl[F[x,y,z],Cartesian[x,y,z]]

Thus, the surface integral is given by

(∫_0^(Pi/2))(∫_0^(2Pi))curlF[S[u,v][[1]],S[u,v][[2]],S[u,v][[3]]] .n[u,v]\[DifferentialD]u\[DifferentialD]v

This answer agrees with the one obtained using the line integral definition.

Example 17.5.  Find the flux of the curl of the vector field F(x,y,z)=⟨ (x^2), (z^2) , y+ (x^2)⟩ across S, where S is the part of the cone (z^2)=(x^2)+(y^2) for which 1≤z≤4 with outward normal orientation.

Solution:  First, we will need the following parametric equations to describe the cone S: x=u cos v , y=u sin v, and z=u, where 0≤v≤2π and 1≤u≤4.

For the cone to have outward orientation, we set n=(Subscript[T, v])×(Subscript[T, u]) (right-hand rule) since (Subscript[T, v]) points in the horizontal direction around the cone and (Subscript[T, u]) points in the direction along the length of the cone.

Clear[F,S,u,v,Tu,Tv,n] 
F[x_,y_,z_]={(x^2)+(y^2), x+(z^2),0}
S[u_,v_]:= { u Cos[v]  , u Sin[v] ,u } 
Tu[u_,v_]:=D[S[u,v],u]
Tv[u_,v_]:=D[S[u,v],v]
n[u_,v_]=Cross[Tv[u,v],Tu[u,v]]

We now compute the flux of curl(F) across S through the following steps.

curlF[x_,y_,z_]=Curl[F[x,y,z],Cartesian[x,y,z]]

Flux=(∫_1^4)(∫_0^(2Pi))curlF[S[u,v][[1]],S[u,v][[2]],S[u,v][[3]]].n[u,v]\[DifferentialD]v\[DifferentialD]u

Exercises

NOTE: In order to perform the curl operation in Mathematica, it is necessary to first load the VectorAnalysis package.  See instructions given at the beginning of this chapter.

In Exercises 1 and 2, find the curl of the given vector field.
1. F(x,y,z)=⟨ ln ((x^2)+(y^2)+(z^2)),x/ z ,(e^x)sin (y z)⟩
2. F(x,y,z)=⟨ -(x/((((x^2)+(y^2)+(z^2)))^(3/2))),-(y/((((x^2)+(y^2)+(z^2)))^(3/2))),-(z/((((x^2)+(y^2)+(z^2)))^(3/2)))⟩

In Exercises 3 and 4, verify Stokes's Theorem for the given vector field F and surface S.
3. F(x,y,z)=⟨ (x^3)e-3x y+(z^3),2 (z^3)-x (z^2)+(y^4),6y+2(z^3)(x^2)⟩ and S is the part of the paraboloid z=(x^2)+(y^2) for which z≤9 and with outward normal orientation.
4. F(x,y,z)=⟨ x y z,x y,x+y+z⟩ and S is the elliptical region in the plane y+z=2 whose boundary is the intersection of the plane with the cylinder (x^2)+(y^2)=1 and with upward normal orientation.

In Exercises 5 and 6, use Stokes's Theorem to compute the flux of the curl of the vector field F across the surface S.
5. F(x,y,z)=⟨ tan(x y z),(e^(y-x z)),sec((y^2)x)⟩ and S is the upper hemisphere of radius 4.
6. F(x,y,z)=⟨ (x^2)z,x (y^2),(z^2)⟩ and S consists of the top and four sides of the cube (excluding the bottom) with vertices at (0,0,0), (1,0,0), (0,1,0), (1,1,0), (0,0,1), (1,0,1), (0,1,1), (1,1,1).

17.3  Divergence Theorem

Students should read Section 17.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

Let F(x,y,z)=⟨(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])⟩ be a vector field.  The divergence of F, denoted by div(F) or ∇·F, is defined by

div(F)=∇·F=((∂(Subscript[F, 1]))/(∂x)) +((∂(Subscript[F, 2]))/(∂ y))+((∂(Subscript[F, 3]))/(∂z))

where ∇= ⟨(∂/(∂x)),(∂/(∂y)),(∂/(∂z))⟩.

The Mathematica command for computing the divergence of a vector field F is Div[F,coordsys], where coordsys is the coordinate system of the vector field.  This is demonstrated in the next example.

The following is another generalization of the Fundamental Theorem of Calculus three dimensions, which relates a triple integral of a solid object involving divergence with a corresponding surface integral along its boundary.

Divergence Theorem: Let W be a region in (ℝ^3) whose boundary ∂W is a piecewise smooth surface, oriented so that the normal vectors to ∂W point outside of W, and F(x,y,z) be a vector field with continuous partial derivatives whose domain contains W.  Then

∫(∫_(∂W))F·dS =∫∫(∫_W)div(F) dV

Refer to your textbook for a detailed discussion and proof of the Divergence Theorem.

Example 17.8.  Find the divergence of the vector field F(x,y,z)=⟨ x sin (yz),(e^(x /y)) z ,(yx^2)⟩.

Solution:

Clear[F1,F2,F3,x,y,z]
F1 =x  Sin[y z]
F2  =(E^(x/y))z
F3 =((x )^2) y
F={F1,F2,F3}

Then the divergence of F is

Div[F,Cartesian[x,y,z]]

NOTE: Again we obtain the same answer for the divergence of F using the explicit formula:

D[F1,x]+D[F2,y]+D[F3,z]

Example 17.9.  Find ∫(∫_S) F·dS, where F(x,y,z)=⟨ x ,(y^2),y+z⟩ and S=∂W is the boundary of the region W contained in the cylinder  (x^2)+(y^2)=4 between the plane z=x and z=8.

Solution: If S is the boundary of the solid W, then W is given by

W={(x,y,z): -2≤x≤2, -(Sqrt[4-(x^2)])≤y≤(Sqrt[4-(x^2)]), x≤z≤8}

Clear[F,divF,x,y,z]
F[x_,y_,z_]={x,(y^2),y+z}
divF=Div[F[x,y,z],Cartesian[x,y,z]]

By the Divergence Theorem, we see that ∫(∫_S) F·dS is given by

(∫_(-2)^2)(∫_(-(Sqrt[4-(x^2)]))^(Sqrt[4-(x^2)]))(∫_x^8)divF\[DifferentialD]z\[DifferentialD]y\[DifferentialD]x

Exercises

NOTE: In order to perform the divergence operation in Mathematica, it is necessary to first load the VectorAnalysis package.  See instructions given at the beginning of this chapter.

In Exercises 1 and 2, find the divergence of the given vector field F.
1. F(x,y,z)=⟨ x y z,(x^2)+(y^2)+(z^2),x y+y z+x z⟩
2. F(x,y,z)=⟨ (e^(x y))cos z,(e^(y z))sin z, (z^2)⟩

In Exercises 3 and 4, verify the Divergence Theorem for the given vector field F and solid region W.
3. F(x,y,z)=⟨(x^2)y,(y^2)z,(z^2)x ⟩ and W={(x,y,z):(x^2)+(y^2)+(z^2)<1} is the unit ball. 
4. F(x,y,z)=⟨(e^x)cos y,(e^x)sin y,x y z ⟩ and W is the region bounded by the paraboloid z=(x^2)+(y^2) and z=4.

In Exercises 5 and 6, use the Divergence Theorem to calculate the flux of the vector field F across the surface S.
5. F(x,y,z)=⟨ x (e^(z )),(y^2),y+z x⟩ and S is tetrahedron bounded by the plane 3x+4y+5z=15 and the coordinate planes in the first octant.
6. F(x,y,z)=⟨ x y z,(x^2)+(y^2)+(z^2),x y+y z+x z ⟩ and S is the unit cube with vertices at (0,0,0), (1,0,0), (0,1,0), (1,1,0), (0,0,1), (1,0,1), (0,1,1), (1,1,1).