Example: Let's define the vector field and test to see if it's conservative by computing . The curl of will be the zero vector field if is conservative.  

>	F := VectorField(<x^2+y^2,2*x*y,y^2-4*x*z>); Curl(F);

 

 

	(1)

 

Since is not zero, we conclude that is not conservative. However, the divergence of is zero. 

>	Divergence(F);

 

(2)

 

Next, let  

>	F := VectorField(<y*z+3*x^2*y^2+2*x,x*z+2*x^3*y+z*cos(y),x*y+sin(y)>); Curl(F); Divergence(F);

 

 

 

	(3)

 

This time . Let's find a potential function for . 

>	f := ScalarPotential(F); 'F(x,y,z)' = Gradient(f,[x,y,z]);

 

 

	(4)

 

Now consider the line integral  where is the upper-half of the circle lying in the xy-plane centered at (2,0,0) with radius 2 oriented counter-clockwise. Let's compute this integral three different ways. First, we'll compute it directly after parametrizing as follows: , , and where . 

>	r := <2*cos(t)+2,2*sin(t),0>;

 

(5)

 

We plug our parametrization into our vector field using the "evalVF" (evaluate vector field command) and then dot the result with the derivative of our parametrization . Then we must integrate from to (half way around the circle). 

>	Int(evalVF(F,r).diff(r,t),t=0..Pi) = int(evalVF(F,r).diff(r,t),t=0..Pi);

 

(6)

 

For a second try at evaluating this integral, remember that our vector field is conservative. So we can use the fundamental theorem of line integrals to find the answer. We all ready have a potential function. We also need to find the start and end points of our curve .  

>	'r(0)' = subs(t=0,r); 'r(Pi)' = subs(t=Pi,r);

 

 

	(7)

 

Thus the start of is (4,0,0) and the end is (0,0,0). 

>	'f(0,0,0)' - 'f(4,0,0)' = subs({x=0,y=0,z=0},f) - subs({x=4,y=0,z=0},f);

 

(8)

 

Now let's use compute the integral using the relation:   .  We need to compute both (the unit tangent for ) and . Recall that . 

>	ds := simplify(Norm(diff(r,t))); T := diff(r,t)/Norm(diff(r,t)); Int(evalVF(F,r).T * ds,t=0..Pi) = int(evalVF(F,r).T * ds,t=0..Pi);

 

 

 

	(9)

 

Example: First, we'll define a vector field and a parametric curve where then we'll plot our curve. 

>	F := VectorField(<3*x^2*y^2*z+y, 2*x^3*y*z+x, x^3*y^2+cos(z)>); r := <cos(t),sin(t),t>; spacecurve(r,t=0..3*Pi/2,scaling=constrained,axes=boxed);

 

Next, we compute the line integral  where is the curve parametrized by . The command "evalVF(F,r)" plus our parameterization into our vector field (it evaluates at ). Then we dot this with the derivative of our parametrization and integrate. 

>	int(evalVF(F,r).diff(r,t),t=0..3*Pi/2);

 

(10)

 

Let's recompute the line integral  by computing . [Of course, we should get the same answer.] First, recall that and is the unit Tangent function. 

>	ds := Norm(diff(r,t)); T := Normalize(diff(r,t)); int((evalVF(F,r).T)*ds,t=0..3*Pi/2);

 

 

 

	(11)

 

Next, let's verify that is conservative and then use the fundamental theorem of line integrals to compute . 

>	Curl(F); f := ScalarPotential(F); startPt := subs(t=0,r);  endPt := subs(t=3*Pi/2,r); simplify(subs({x=0,y=-1,z=3*Pi/2},f) - subs({x=1,y=0,z=0},f));

 

 

 

 

 

	(12)

 

Example: Find the centroid of : the upper-half of the ellipse .  We know that the centroid is given by where ,  , and . First, we will parametrize and then plot . 

>	x := t -> 3*cos(t): y := t -> 2*sin(t): r := t -> <x(t),y(t)>: 'r(t)' = r(t); plot([x(t),y(t),t=0..Pi],scaling=constrained,color=blue,thickness=3);

 

Now let's find the arc length "" of and then find the moments about the y and x axes and compute the centroid. 

>	ds := Norm(diff(r(t),t)); m := int(1*ds,t=0..Pi); M_y := int(x(t)*ds,t=0..Pi); M_x := int(y(t)*ds,t=0..Pi); Centroid := [M_y/m, M_x/m]; evalf(Centroid);

 

 

 

 

 

 

	(13)

 

>	# Let's undefine x,y,z for later examples: x := 'x': y := 'y': z := 'z':

 

Example: Let C be the upper-half of the circle which lies in the plane oriented counter-clockwise (beginning at (3,0,0) and ending at (-3,0,0)). Let . First, let's compute the line integral by plugging in the a parameterization for C and integrating. 

>	F := VectorField(<y*z-3*x^2+x^z*z/x,x*z-z*exp(y*z)-10*y,x*y-y*exp(y*z)+x^z*ln(x)>); r := <3*cos(t),3*sin(t),0>; int(evalVF(F,r).diff(r,t),t=0..Pi);

 

 

 

	(14)

 

Now let's recompute the line integral using the fundamental theorem of line integrals -- that is -- first, we find a potential function for and then we plug in the endpoints of C. [I have included a "curl" command which verifies that is indeed conservative.] 

>	Curl(F); f := ScalarPotential(F); subs(x=-3,y=0,z=0,f) - subs(x=3,y=0,z=0,f);

 

 

 

	(15)

 

Example: Let's compute the line integral where  and C is the streched helix defined by and .   

>	P := y*z+y^2+exp(x^2); Q := x*z+2*x*y+cos(z)+1; R := x*y-y*sin(z)+sin(z)/z; F := VectorField([P,Q,R]); r := t -> <2*cos(t),3*sin(t),t^2>: 'r(t)' = r(t); spacecurve(r(t),t=0..4*Pi,axes=boxed,color=blue,thickness=3,orientation=[60,70]);

 

Let's see if is conservative. If so, its curl should be zero...which it is. 

>

Curl(F);

 

(16)

 

Next, let's compute the integral  using the definition of this type of line integral and the given parametrization. We'll use "evalVF" (evaluate vector field) to plug the parametrization into the vector field. Also, we'll use "evalf" to find a decimal approximation since Maple won't compute our integral exactly. 

>	evalf(int(evalVF(F,r(t)).diff(r(t),t),t=0..4*Pi));

 

(17)

 

Now since our vector field is conservative, we can compute the integral  using the fundamental theorem of line integrals. Let's do so. First, we'll compute a potential function. 

>	f := ScalarPotential(F);

 

(18)

 

I'm glad we didn't try to do that one by hand. The potential function involves the error function and the sine integral functions (I didn't learn about either of those in pre-calculus). Next, we need the end points of our curve. 

>	r(0); r(4*Pi);

 

 

	(19)

 

So the curve starts at and ends at . Now plugging in the end points and computing their difference will yield the answer. 

>	subs({x=2,y=0,z=16*Pi^2},f) - subs({x=2,y=0,z=0},f);

 

(20)

 

The answers should match... 

>	evalf(subs({x=2,y=0,z=16*Pi^2},f) - subs({x=2,y=0,z=0},f));

 

(21)

 

And so they do. Considering that , we have that the exact value of our line integral is . It is interesting to note that we can get an "exact" answer (even if it involves an unfamilar function) by using the fundamental theorem of line integrals while computing the integral directly yields an integral that Maple cannot handle exactly. This is kind of like the multiple integrals that we ran into before where Maple could handle the integral in cylindrical or spherical coordinates, but not in rectangular coordinates. 

Example: For each of the following vector fields, let's determine if it is either the curl or gradient of something.  

If you look at the vector fields defined below, they are built from functions which have continuous derivatives (of all orders) everywhere. For "nice" vector fields like these the following theorems hold: if and only if  for some scalar valued function . In other words, vector fields with potential functions (conservative vector fields) have zero curl ("Curl( Gradient( f ) )=0") and vector fields with zero curl are conservative. [We should note that the second half of this requires the vector field to be "nice" on a simply connected domain. This is not a problem for us since our vector fields are "nice" everywhere.] The second theorem (which isn't important to us -- in fact -- it's not even discussed in our text) says that if and only if for some vector field . In other words, vector fields which are the result of curl-ing another vector field have zero divergence ("Divergence( Curl( F ) )=0") and vector fields with zero divergence are the curl of some vector field. 

In Maple, if a vector field is the gradient of something, we can "un-gradient" it with the function "ScalarPotential". If the vector field is the curl of something, we can "un-curl" it with the function "VectorPotential". 

>	F := VectorField(<2*x+y^2*exp(x*y^2), 2*x*y*exp(x*y^2)+2*sin(y+2*z)*cos(y+2*z), 4*sin(y+2*z)*cos(y+2*z)+exp(-z)-z*exp(-z)>); 'curl' = Curl(F); 'divergence' = Divergence(F); f := ScalarPotential(F); 'F' = Gradient(f,[x,y,z]);

 

 

 

 

 

	(22)

 

F has no vector potential since div(F) is not 0. However, it did have a (scalar) potential function since curl(F)=0. 

>	G := VectorField(<x+y+z,x^3*y,-x^3*z-z>); 'curl' = Curl(G); 'divergence' = Divergence(G); M := VectorPotential(G); 'G' = Curl(M);

 

 

 

 

 

	(23)

 

has no scalar potential function since . However, it does have a vector potential (so ). 

>	H := VectorField(<x*y,y*z,x*z>); 'curl' = Curl(H); 'divergence' = Divergence(H);

 

 

 

	(24)

 

has no vector or scalar potential function since and . 

You might ask, "What happens if we ask Maple to find a vector or scalar potential when these do not exist?" The answer is "nothing" (no output at all). 

>	VectorPotential(H); ScalarPotential(H);

 

Example: [Green's Theorem] Let be the curve which consists of the parabola from to and then the line segment from to. Let's plot and define the functions:  and . 

>	display(plot(x^2,x=0..1),plot(x,x=0..1));

 

>	P := (x,y) -> y^2*sin(x): 'P(x,y)' = P(x,y); Q := (x,y) -> x^2*sin(y): 'Q(x,y)' = Q(x,y);

 

 

	(25)

 

Let's compute the line integral . To do this we need to parameterize both pieces of our curve, plug everything in for each curve, integrate along both pieces, and then put the results together. Notice that for parametrizes the parabolic piece and for parametrizes the line segment. 

>	x_1 := t; y_1 := t^2; x_2 := 1-t; y_2 := 1-t; A := int(P(x_1,y_1)*diff(x_1,t) + Q(x_1,y_1)*diff(y_1,t), t=0..1); B := int(P(x_2,y_2)*diff(x_2,t) + Q(x_2,y_2)*diff(y_2,t), t=0..1); lineInt := simplify(A+B);

 

>

 

 

 

 

 

 

 

	(26)

 

Now let's use Green's Theorem to turn the line integral above into an equivalent double integral (and verify that Green's Theorem holds): . 

Notice that our curve is closed and oriented in a counter-clockwise direction. It is easy to see that and in our region (the region bounded by ). 

>	doubleInt := int(int(diff(Q(x,y),x)-diff(P(x,y),y),y=x^2..x),x=0..1); evalb(doubleInt = lineInt);

 

 

	(27)

 

Example: We will find the area enclosed by the "bow-tie" curve C which is parametrized by and where . Let's find this area using one of the line integral formulas (which follow from Green's theorem):  . Note: These area formulas assume that the curve has been traversed in a counter-clockwise direction (note the animation below). If we're not careful to pay attention to the orientation of our curve we'll get a wrong answer. 

>	plot([3*sin(t),2*sin(2*t),t=0..2*Pi],scaling=constrained);

 

After executing the next command, click on the the picture below. A toolbar with animation controls will appear. 

>	animate([3*sin(c*t),2*sin(2*c*t), t=0..2*Pi], c=0..1, view=[-3..3,-2..2], numpoints=150, frames=64, scaling=constrained);

 

Keep in mind that the line integral only computes area if the curve goes around that region in a counter-clockwise direction. If the curve goes clockwise, then the line integral(s) will compute the negative of the area. That's why with this curve if you integrate from 0 to 2π the answer is 0 since half of the "bow-tie" goes around clockwise and half goes counter-clockwise.  

Therefore, to compute the area we should integrate from π to 2π (where the curve is travelling in the counter-clockwise direction and then (using symmetry) double the result to find out area. 

Here are a few different ways we could arrive at the our answer of "16"... [the first three integrals just find the area of the left half (counter-clockwise orientated half) of the bow-tie and double the answer. The last line computes the integral over the right half and negates (since the curve is oriented clockwise on this half) and then adds in the left half. 

>	x := 3*sin(t); y := 2*sin(2*t); dx := diff(x,t); dy := diff(y,t); "Area" = 2*int(x*dy, t=Pi..2*Pi); "Area" = 2*int(-y*dx, t=Pi..2*Pi); "Area" = 2*int((1/2)*(x*dy-y*dx), t=Pi..2*Pi); "Area" = -int(x*dy, t=0..Pi) + int(x*dy, t=Pi..2*Pi);

 

 

 

 

 

 

 

 

	(28)

 

Now let's compute the arc length of this curve. That is find:  . Note:   

>	ds := sqrt(dx^2+dy^2); "Arc Length" = evalf(int(1*ds,t=0..2*Pi));

 

 

	(29)

 

Note: The "evalf" command tells Maple to approximate the answer. The integral is so "hard" that Maple can't get an exact answer. 

Example: Let's find the area inside the following curve using a line integral area formula (which follow from Green's theorem): etc.  

>

x := t -> 1/10*cos(t)+cos(t/10): y := t -> 1/10*sin(t)-sin(t/10): 'r(t)' = <x(t),y(t)>; plot([x(t),y(t),t=0..20*Pi],scaling=constrained,color=green,thickness=3); plotA := plot([x(t),y(t),t=0..2*Pi],color=blue,thickness=5): plotB := plot([x(t),y(t),t=0..20*Pi],color=green,thickness=1): display({plotA,plotB},scaling=constrained);

 

The second plot shows the whole curve (in green) along with just the portion where (in blue). We can see that our curve is being traced out in the clockwise direction. I'll use to compute the area... 

>	int(x(t)*diff(y(t),t),t=0..20*Pi);

 

(30)

 

Notice that we got a negative answer. This is due to the fact that the curve is oriented in a clockwise direction (the "wrong" direction). But no worries, merely flipping the sign will fix the answer. Thus the area inside the curve is . 

>	# Let's undefine x,y,z for later examples: x := 'x': y := 'y': z := 'z':

 

A bit about Orientations 

We deal with two different types of surfaces in class. (1) Level surfaces defined by where is some fixed constant. (2) Parametrized surfaces defined by where belong to some region . Mostly, we deal with graphs of functions . These fit into both (1) and (2). Notice that such a graph is a level curve and it can be parameterized by .  

An orientation for a surface is a smooth choice of unit normal vectors - that is an orientation is a vector field on a surface such that at each point the vector field describes a (unit length) normal vector for the tangent plane of the surface. 

We know that gives us normal vectors for the tangent planes of a level surface . So (as long as  ≠ 0 ) we get the orientations: and for our surface . 

If our surface is parametrized by we can find tangent vectors for our surface by computing and . So gives normal vectors for the planes tangent to our surface (as long as isn't zero). So (as long as ) we get the orientations: n = and n = . 

Example: Consider the surface defined by where . We can view as the level surface 0 where . This gives us the following orientations: 

>	F := x^2+y^2-z; gradF := Gradient(F,[x,y,z]); nDown := Normalize(gradF); nUp := -Normalize(gradF);

 

 

 

 

	(31)

 

We can distinguish between the downward and upward orientations by looking at the "z" component of the orientation. Notice that this component in nDown is negative wheras the corresponding component in nUp is positive. 

Next, we ca also view as a parametrized surface using the parametrization: ,  ,  where . This gives us the following orientations: 

>	r := <u,v,u^2+v^2>; rU := diff(r,u); rV := diff(r,v); nUp := Normalize(rU &x rV); nDown := -Normalize(rU &x rV);

 

 

 

 

 

	(32)

 

Let's graph along with it's downward orientation. [You may ask, "What are those commands and all those complicated formulas for?" Well, is the sequence command. This allows us to plug in a sequence of values into a formula. The whole idea of those commands is to list off a bunch of sample normal vectors for our surface without having to list them individually.] 

>

plotS := plot3d(r,u=-2..2,v=-sqrt(4-u^2)..sqrt(4-u^2),numpoints=1000): normalPlot := seq(seq(arrow(subs({u=a,v=b/((4-a^2)+0.01)*sqrt(4-a^2)},r),                            subs({u=a,v=b/((4-a^2)+0.01)*sqrt(4-a^2)},nDown),                            shape=arrow,color=black),b=-(4-a^2)..(4-a^2)),a=-2..2): display({plotS,normalPlot},scaling=constrained,viewpoint=circleleft);

 

Example: Consider the surface where . We could plot this surface using "implicitplot3d" or by first parameterizing the surface and then using "plot3d". Let's do both.  

>	implicitplot3d(sqrt(x^2+y^2)-z, x=-3..3, y=-3..3, z=1..3, numpoints=5000, scaling=constrained, axes=boxed);

 

Noticing the form of the equation, we should expect that cylindrical coordnates should yield a nice parameterization. In cylindrical coordinates, . So using and as parameters, we see that and (there's no reason to restrict θ). Therefore, we get: , and where and . [We get a nicer plot!] 

>	# undefine r r := 'r'; plot3d([r*cos(theta),r*sin(theta),r],r=1..3,theta=0..2*Pi,scaling=constrained,axes=boxed);

 

Next, let's find formulas for both orientations of this surface using the gradient operator and then find orientations using our parameterization. [We can distinguish between upward and downward orientations by considering the sign on the z-component of our orientation.] 

>	gradF := Gradient(sqrt(x^2+y^2)-z); orientUp := -Normalize(gradF); orientDown := Normalize(gradF);

 

 

 

	(33)

 

Now orientations via the parameterization... 

>	X := (r,theta) -> <r*cos(theta),r*sin(theta),r>: 'X(r,theta)' = X(r,theta); orientUp := Normalize(diff(X(r,theta),r) &x diff(X(r,theta),theta)) assuming r::positive; orientDown := -Normalize(diff(X(r,theta),r) &x diff(X(r,theta),theta)) assuming r::positive;

 

 

 

	(34)

 

Finally, let's plot this surface (using 's parametric plot) along with a few sample normal vectors (using our downward orientation). 

>	plotS := plot3d(X(r,theta),r=1..3,theta=0..2*Pi): normalPlot := plotNormals(X,[1..3,0..2*Pi],[5,15],true): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Example: Let's plot an ellipsoid with its outward orientation. First, notice that we can parameterize an ellipsoid using "modified spherical coordinates": 

>	x^2/a^2+y^2/b^2+z^2/c^2=1; xE := a*rho*cos(theta)*sin(phi); yE := b*rho*sin(theta)*sin(phi); zE := c*rho*cos(phi); simplify(xE^2/a^2+yE^2/b^2+zE^2/c^2=1);

 

 

 

 

 

	(35)

 

In these modified spherical coordinates, we have . This gives the following parametrization for the ellipsoid (allowing and ). Let's also specify a specific ellipsoid (with ). Then we'll compute its outward orientation. 

>	r := (theta,phi) -> <a*cos(theta)*sin(phi), b*sin(theta)*sin(phi), c*cos(phi)>: 'r(theta,phi)' = r(theta,phi); a := 1;  b := 2;  c := 3; nOut := simplify(-Normalize(diff(r(theta,phi),theta) &x diff(r(theta,phi),phi)));

 

 

 

 

 

	(36)

 

Now let's plot our ellipsoid along with a grid of normal vectors. Since the domain for our parameterized surface is a rectangle (in θ,φ space), we can use the plotNormals function I defined at the beginning of the worksheet. 

>	plotS := plot3d(r(theta,phi),theta=0..2*Pi,phi=0..Pi): normalPlot := plotNormals(r,[0..2*Pi,0..Pi],[10,10],true): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Example: Let's plot the part of the unit sphere in the first octant along with some normal vectors. Then we'll compute its surface area and centroid. 

>	r := (u,v) -> <cos(u)*sin(v),sin(u)*sin(v),cos(v)>: 'r(u,v)'=r(u,v);     plotS := plot3d(r(u,v),u=0..Pi/2,v=0..Pi/2,scaling=constrained): normalPlot := plotNormals(r,[0..Pi/2,0..Pi/2],[7,5]): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Adding the optional argument "true" to the end of the "plotNormals" command will flip the orientation and show outward pointing normals. 

>	plotS := plot3d(r(u,v),u=0..Pi/2,v=0..Pi/2,scaling=constrained): normalPlot := plotNormals(r,[0..Pi/2,0..Pi/2],[7,5],true): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Now let's compute the surface area and centroid of this surface. 

>	dS := simplify(Norm(diff(r(u,v),u) &x diff(r(u,v),v))); m := int(int(1*dS,u=0..Pi/2),v=0..Pi/2); Myz := int(int(cos(u)*sin(v)*dS,u=0..Pi/2),v=0..Pi/2); Mxz := int(int(sin(u)*sin(v)*dS,u=0..Pi/2),v=0..Pi/2); Mxy := int(int(cos(v)*dS,u=0..Pi/2),v=0..Pi/2); Centroid = [Myz/m,Mxz/m,Mxy/m];

 

 

 

 

 

 

	(37)

 

Example: Lets find the centroid of the part of the surface which lies above the -plane. First, we'll plot our surface. The bounds for and can be found by noticing that their extreme values are determined by where this surface intersects the -plane (i.e. ). Thus  solving for we get . This is an ellipse. The extreme values for can be found by setting these equations to zero. We get . 

>	plot3d(1-x^2/9-y^2/16,y=-4*sqrt(1-x^2/9)..4*sqrt(1-x^2/9),x=-3..3,scaling=constrained,axes=boxed);

 

We could parametrize our surface using: , , and where and . But this is quite messy. Instead let's use modified polar coordinates (sort of "elliptic coordinates"). Keeping in mind that we want to clear the 9 and 16 in the denominators, we should scale up by 3 and by 4. This gives us , , and .  Again determines extreme values, so tells us that (and ).  As to not overlap with the name "r" which we tend to use to name parametrizations, I will use parameters and . 

>	x := (u,v) -> 3*u*cos(v): y := (u,v) -> 4*u*sin(v): z := (u,v) -> 1-u^2: r := (u,v) -> <x(u,v),y(u,v),z(u,v)>: 'r(u,v)' = r(u,v);

 

(38)

 

Let's plot our surface again. This time we'll use our "nice" parametrization (notice the change in the grid lines on our surface). 

>	plot3d(r(u,v),u=0..1,v=0..2*Pi,scaling=constrained,axes=boxed);

 

To compute the centroid of our surface we need to compute several surface integrals. All such integrals involve the surface area element "dS". By definition: .   

>	dS := Norm(diff(r(u,v),u) &x diff(r(u,v),v));

 

(39)

 

Now we need to compute the mass (i.e. the surface area) and 3 moments:,  ,   ,   and   . 

>	m := evalf(int(int(dS,u=0..1),v=0..2*Pi)); M_yz := evalf(int(int(x(u,v)*dS,u=0..1),v=0..2*Pi)); M_xz := evalf(int(int(y(u,v)*dS,u=0..1),v=0..2*Pi)); M_xy := evalf(int(int(z(u,v)*dS,u=0..1),v=0..2*Pi));

 

 

 

 

	(40)

 

I chose to use "evalf" since the exact values of these integrals give Maple trouble. Notice that and are nearly 0. In fact, we can tell by symmetry that the and coordinates of the centroid should be exactly 0. So we'll replace these values with 0 and chalk up the difference to round-off error. 

>	Centroid := [0,0,M_xy/m];

 

(41)

 

>	# Let's undefine x,y,z for later examples: x := 'x': y := 'y': z := 'z':

 

Example: Consider the torus parameterized by where and . Let's graph the torus and find its surface area. 

>	r := (u,v) -> <(3+cos(v))*cos(u),(3+cos(v))*sin(u),sin(v)>: 'r(u,v)' = r(u,v); plot3d(r(u,v),u=0..2*Pi,v=0..2*Pi,scaling=constrained,viewpoint=circleleft); dS := simplify(Norm(diff(r(u,v),u) &x diff(r(u,v),v))); Int(Int(1*dS,u=0..2*Pi),v=0..2*Pi) = int(int(1*dS,u=0..2*Pi),v=0..2*Pi);

 

 

 

 

	(42)

 

In general, a torus is built by rotating a circle of radius "" around a circle of radius "". The resulting torus has surface area . In our example, and so the surface area is . 

Example: Consider the following parameterized surface (it's called a Mobius strip)... 

>	X := (u,v) -> [-cos(u/4)+4*cos(u) + cos(u)*cos(u/2)*v, 4*sin(u) + sin(u)*cos(u/2)*v, sin(u/2)*v]: 'X(u,v)'=X(u,v); XDomain := [0..2*Pi,0..1];

 

 

	(43)

 

The following plot illustrates the fact that a Mobius strip is not orientable (notice how the normal vectors don't match at the "seem"). 

>	mobiusPlot := plot3d(X(u,v),u=XDomain[1],v=XDomain[2]): normalPlot := plotNormals(X,XDomain,[20,2]): display(mobiusPlot,normalPlot,scaling=constrained,orientation=[0,75],viewpoint=circleleft);

 

Although we cannot orient the mobius strip, we can still compute its surface area. Remember that "dS" does not care about orientation! Let's compute the surface area of this surface by first computing "" and then integrating. 

WARNING: Maple can't handle the double integral which computes this surface area (directly). We must use an "evalf" command to get a decimal approximation. Also,  it's best to choose the order of integration: "du dv" not "dv du" -- Maple doesn't seem to be able to handle the order "dv du" even when approximating! 

>	Xvec := convert(X(u,v),Vector): dS := simplify(Norm(diff(Xvec,u) &x diff(Xvec,v))); evalf(int(int(1*dS,u=0..2*Pi),v=0..1));

 

 

	(44)

 

Example: Consider the representation of the Mobius strip given by where and . 

>	r := (u,v) -> <2*cos(u)+v*cos(u/2),2*sin(u)+v*cos(u/2),v*sin(u/2)>: 'r(u,v)' = r(u,v);     plotS := plot3d(r(u,v),u=0..2*Pi,v=-1/2..1/2): normalPlot := plotNormals(r,[0..2*Pi,-1/2..1/2],[20,2]): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

We can see from the plot above that the Mobius strip is not orientable. However, we can still find its surface area and centriod. 

>

dS := simplify(Norm(diff(r(u,v),u) &x diff(r(u,v),v))); m := int(int(1*dS,u=0..2*Pi),v=-0.5..0.5); Myz := int(int((2*cos(u)+v*cos(u/2))*dS,u=0..2*Pi),v=-0.5..0.5); Mxz := int(int((2*sin(u)+v*cos(u/2))*dS,u=0..2*Pi),v=-0.5..0.5); Mxy := int(int((v*sin(u/2))*dS,u=0..2*Pi),v=-0.5..0.5); Centriod = [Myz/m,Mxz/m,Mxy/m];

 

 

 

 

 

 

	(45)

 

Example: The following parametrization is part of a torus (a doughnut shaped surface). 

>	r := (u,v) -> <(3+cos(v))*cos(u),(3+cos(v))*sin(u),sin(v)>: 'r(u,v)' = r(u,v);     plotS := plot3d(r(u,v),u=0..Pi,v=0..Pi,scaling=constrained): normalPlot := plotNormals(r,[0..Pi,0..Pi],[15,5]): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Let's compute the surface area and centroid for this surface. 

>	dS := simplify(Norm(diff(r(u,v),u) &x diff(r(u,v),v))); m := int(int(1*dS,u=0..Pi),v=0..Pi); Myz := int(int((3+cos(v))*cos(u)*dS,u=0..Pi),v=0..Pi); Mxz := int(int((3+cos(v))*sin(u)*dS,u=0..Pi),v=0..Pi); Mxy := int(int(sin(v)*dS,u=0..Pi),v=0..Pi); Centroid = [Myz/m,Mxz/m,Mxy/m];

 

 

 

 

 

 

	(46)

 

Example: [Stokes' Theorem] Let be the part of the paraboloid which lies above the disk . Also, let be oriented "upward" (the component of the orientation should be positive). Recall that (the boundary of our surface) inherits the counter-clockwise orientation since our surface is oriented upwards. Also, let . 

Let's compute the flux integral by plugging in a parameterization for the surface. Notice that this surface is most easily described in terms of cylindrical coordinates: and . So we get the following parameterization for our surface: and  where and . We'll begin by plotting our surface (with some of it's normal vectors). 

>	# undefine "r" (since we're using it as a parameter instead of a parameterization name.) r := 'r':

 

>	X := (r,theta) -> [r*cos(theta),r*sin(theta),4-r^2]: 'X(r,theta)' = X(r,theta);     plotS := plot3d(X(r,theta),r=0..2,theta=0..2*Pi): normalPlot := plotNormals(X,[0..2,0..2*Pi],[10,10]): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Let's define our vector field and parameterization (both as vectors). 

>	F := VectorField(<y*z,x*z^3+x*exp(z),y^3*z>); X := <r*cos(theta),r*sin(theta),4-r^2>;

 

 

	(47)

 

Let's compute to see if this yields upward or downward pointing normal vectors. 

>	simplify(diff(X,r) &x diff(X,theta));

 

(48)

 

Since the z-component is and , these are upward pointing. Now we are ready to plug everything in and compute our flux integral. 

>	int(int(evalVF(Curl(F),X).(diff(X,r) &x diff(X,theta)),r=0..2),theta=0..2*Pi);

 

(49)

 

Now let's verify Stokes' theorem by computing the line integral (we should get again).  The boundary of the paraboloid occurs when (i.e. ). This means . Keeping in mind the induced orientation is "counter-clockwise" we get the parameterization: where . 

>	X := <2*cos(t),2*sin(t),0>; int(evalVF(F,X).diff(X,t),t=0..2*Pi);

 

 

	(50)

 

So we found just as Stokes predicts. 

Example: [Divergence Theorem] Let be the spherical region and be the unit sphere (with the induced "outward" orientation). Let . First, let's parameterize our boundary surface . Since this is the unit sphere, spherical coordinates are a natural choice. We have where and . Let's plot this surface with a few outward pointing normal vectors. 

>	X := (theta,phi) -> [cos(theta)*sin(phi),sin(theta)*sin(phi),cos(phi)]: 'X(theta,phi)' = X(theta,phi); plotS := plot3d(X(theta,phi),phi=0..Pi,theta=0..2*Pi): normalPlot := plotNormals(X,[0..2*Pi,0..Pi],[10,10],true): display(plotS,normalPlot,scaling=constrained,viewpoint=circleleft);

 

Next, let's define our vector field and parameterization (as vectors) and then verify that is a formula for outward pointing normal vectors. 

>	F := VectorField(<x^3+y,y^3+x*z,z^3+5*x^2*y^3>); X := <cos(theta)*sin(phi),sin(theta)*sin(phi),cos(phi)>; simplify(diff(X,phi) &x diff(X,theta));

 

 

 

	(51)

 

If we stare at the result of computing    for a while, we should be able to convince ourselves that these normals are outward pointing. Finally, let's go ahead and compute the flux integral . 

>	int(int(evalVF(F,X).(diff(X,phi) &x diff(X,theta)),theta=0..2*Pi),phi=0..Pi);

 

(52)

 

Now verify the divergence theorem by computing the triple integral . Of course, we should compute this triple integral using spherical coordinates (since is the interior of the unit sphere). 

>	Divergence(F);

 

(53)

 

Therefore, .  

>	int(int(int(3*rho^2*rho^2*sin(phi),rho=0..1),theta=0..2*Pi),phi=0..Pi);

 

(54)

 

So we get just as the divergence theorem predicts.