<Text-field style="_cstyle290" layout="Heading 1"><Font bold="true">Visualizing and Interpreting the Divergence Theorem</Font></Text-field>

<Text-field style="_cstyle290" layout="Heading 1"><Font bold="true">Visualizing and Interpreting the Divergence Theorem</Font></Text-field>

Note: You may notice differences between this Maple worksheet and the equivalent Mathematica notebook. These differences were introduced to preserve the content of these modules and were necessary because of major functional differences between Maple and Mathematica.

<Text-field style="Heading 1" layout="Heading 1">Introduction</Text-field>

OBJECTIVE: Observe that surface integrals are difficult to evaluate, even using Maple. See that using parametrizations to evaluate flux surface integrals and applying the Divergence Theorem can help with integral evaluations.

You will see that surface integrals are still difficult to set up and to evaluate, even with Maple, but parametrizations can assist in evaluating these difficult flux surface integrals. You will also verify the Divergence Theorem, just as you have done with Green's Theorem. This module is an example. However, those familiar with Maple can adjust the code to solve other problems.

<Text-field style="_cstyle256" layout="Heading 1"><Font size="18">Technology Guidelines</Font></Text-field>

NOTE: If you have just finished a worksheet, restart Maple before executing a new worksheet. TO OPEN SECTIONS, Click on the Arrow sign at the left hand side of the screen or select Expand All Sections from the View drop down menu.

TO STOP AN EXECUTION Click on STOP button from the toolbar.

ORDER OF EXECUTION Execute commands in the order given. Do not skip any Maple Input lines within a given worksheet

Alternatively, you can execute the entire worksheet by selecting the Execute Worksheet command from the Edit drop down menu.

SAVING WORKSHEETS.

You can save anytime to any directory you choose, and it is wise to save often. EXPERIENCING MAJOR PROBLEMS Save if appropriate, and then shut down Maple and start it up.

<Text-field style="Heading 1" layout="Heading 1">Part I: Visualizing the Problem</Text-field>

We select the object formed by a hemisphere of radius 1, topped with a cylinder of radius 1 and height 1, with a circular top on top on the cylinder. The object is similar to the one to which the his is the method of moments was applied in the previous chapter. The force we choose has components LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNiotSSNtaUdGJDYjUSFGJy1GIzYoLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y/LyUpc3RyZXRjaHlHRj8vJSpzeW1tZXRyaWNHRj8vJShsYXJnZW9wR0Y/LyUubW92YWJsZWxpbWl0c0dGPy8lJ2FjY2VudEdGPy8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRk4tRiM2JC1JJW1zdXBHRiQ2JS1GMTYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjtRJ2l0YWxpY0YnLUkjbW5HRiQ2JFEiMkYnRjovJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjotRjc2LVEoJm1pbnVzO0YnRjpGPUZARkJGREZGRkhGSkZMRk8tRiM2KC1GaW42JFEiNEYnRjotRjc2LVExJkludmlzaWJsZVRpbWVzO0YnRjpGPUZARkJGREZGRkhGSi9GTVEmMC4wZW1GJy9GUEZbcEZWRmdvLUYxNiVRInlGJ0ZZRmZuRjpGMEY6LUY3Ni1RIixGJ0Y6Rj0vRkFGZW5GQkZERkZGSEZKRmpvL0ZQUSwwLjMzMzMzMzNlbUYnLUYjNiZGNi1GIzYoLUZpbjYkUSI2RidGOkZnb0ZdcEZnby1GMTYlUSJ6RidGWUZmbkY6RjBGOkZgcC1GIzYmLUZpbjYkUSMxMkYnRjpGZ29GXXFGOkYwRjpGOi8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Y6.

restart;

force:=[-x^2-4*x*y, -6*y*z, 12*z];

Before we apply the Divergence Theorem, we plot the figure. We need to load in a graphing packages to show our plots.

with(plots):

forceplot:=fieldplot3d(force, x=-1.5..1.5, y=-1.5..1.5, z=-1.5...1.5, axes=BOXED, grid=[6,6,6], arrows=THICK): hemisplot:=sphereplot(1, r=0..2*Pi, theta=0..Pi, labels=["x","y","z"], view = [-1..1,-1..1,-1..0], axes=boxed): pcyl:=plottools[cylinder]([0,0,1],1,.01): #cylinder:=cylinderplot(1, r=0..2*Pi, theta=0..2*Pi, view=[-1..1,-1..1,-1..1]): topplot := plot3d([cos(t),sin(t),u],t=-2*Pi..Pi, u=0..1): display({forceplot, hemisplot, pcyl, topplot}, view=[-1.5..1.5, -1.5..1.5, -1.5..1.5]);

<Text-field style="Heading 1" layout="Heading 1">Part II: Evaluating the Divergence Integral</Text-field>

The divergence integral over the volume is easy to set up and calculate. First, we begin by loading in a package that performs dot products and cross products.

with(linalg):

unassign('x','y','z'); divf:=diff(force[1], x)+diff(force[2],y)+diff(force[3],z): print(`The divergence of this force is`, divf); int(int(int(divf, z=-sqrt(1-x^2-y^2)..1), y=-sqrt(1-x^2)..sqrt(1-x^2)), x=-1..1): totalflux:=evalf(%): print(`The value of the divergence integral is`,totalflux);

<Text-field style="Heading 1" layout="Heading 1">Part III: Finding Surface Integrals</Text-field>

The surface integrals evaluating the flux across each surface require considerable set up. As you recall, parametrizations are very convenient in determining line integrals. Likewise, parametrizations can facilitate the computation of surface integrals. We use spherical coordinates to write the equation of the sphere parametrically.

<Text-field style="Heading 1" layout="Heading 1">Hemisphere</Text-field>

unassign('x','y','z','r','theta','phi'); x:=sin(theta)*sin(phi): y:=cos(theta)*sin(phi): z:=cos(phi): r:=[x,y,z]: rtheta:=diff(r,theta): rphi:=diff(r,phi): cr:=linalg[crossprod](rtheta, rphi): integrand:=12*cos(phi)^2*sin(phi)-6*cos(theta)^2*cos(phi)*sin(phi)^3-sin(theta)^2*(4*cos(theta)+sin(theta))*sin(phi)^4: print(`The integrand for the hemispherical surface integral is`, integrand); hemisflux:=int(int(integrand, phi=(Pi/2)..Pi), theta=0..2*Pi): print(`The value of the hemispherical surface integral is`, hemisflux,`that is`,evalf(hemisflux));

The equations of the cylinder and the top can be written using polar coordinates, with z and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvRjNRJ2l0YWxpY0Yn the parameters for the cylinder.

<Text-field style="_cstyle266" layout="Heading 1"><Font size="18">Cylinder</Font></Text-field>

unassign('x','y','z','r'); x:=cos(theta): y:=sin(theta): z:=z: r:=[x,y,z]: rtheta:=diff(r,theta): rz:=diff(r,z): cr:=linalg[crossprod](rtheta, rz): integrand:=-cos(theta)^3-4*cos(theta)^2*sin(theta)-6*z*sin(theta)^2: print(`The integrand for the cylindrical surface integral is`, integrand); cylflux:=int(int(integrand, z=0..1), theta=0..2*Pi): print(`The value of the cylindrical surface integral is`, cylflux, `that is`, evalf(cylflux));

<Text-field style="Heading 1" layout="Heading 1">Top</Text-field>

Polar coordinates r and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg== will be the parameters for the top. Recall that for the top, the normal is always in the increasing z direction. Also the delta area form of the area must contain the polar coordinate r, which varies from 0 to 1.

unassign('x','y','z','r'); x:=cos(theta): y:=sin(theta): z:=1: r:=[x,y,z]: norm1:=[0,0,1]: integrand:=linalg[dotprod](force,norm1): print(`The integrand for the top surface integral is`, integrand); topflux:=int(int(R*integrand, R=0..1),theta=0..2*Pi): print(`The value of the surface integral for the top is`, topflux, `that is` ,evalf(topflux));

<Text-field style="Heading 1" layout="Heading 1">Part IV: Verifying the Divergence Theorem</Text-field>

Compare your answers to verify the Divergence Theorem.

totalsurfaceflux := hemisflux + cylflux + topflux; print(`The total of the surface integrals is `, totalsurfaceflux, `that is`,evalf(totalsurfaceflux)); if (evalf(totalsurfaceflux)<=(totalflux+0.001) and evalf(totalsurfaceflux)>=totalflux-.001) then print(`The Divergence Theorem is validated.`) else print(`there's a problem`) fi;

<Text-field style="Heading 1" layout="Heading 1">You Try It: Part IV</Text-field>

<Text-field style="_cstyle268" layout="Heading 1"><Font size="18">Visualize and Evaluate Other Divergence Integrals</Font></Text-field>

You can change the force function to define a different force function and then to visualize it.

unassign('x','y','z', 't', 'r'); force:=[cos(10*z), sin(10*z), x*y]; forceplot:=fieldplot3d(force, x=-1.5..1.5, y=-1.5..1.5, z=-1.5...1.5, axes=BOXED, grid=[6,6,6], arrows=THICK): hemisplot:=sphereplot(1, r=0..2*Pi, theta=0..Pi, labels=["x","y","z"], view = [-1..1,-1..1,-1..0], axes=boxed): pcyl:=plottools[cylinder]([0,0,1],1,.01): #cylinder:=cylinderplot(1, r=0..2*Pi, theta=0..2*Pi, view=[-1..1,-1..1,-1..1]): topplot := plot3d([cos(t),sin(t),u],t=-2*Pi..Pi, u=0..1): display({forceplot, hemisplot, pcyl, topplot}, view=[-1.5..1.5, -1.5..1.5, -1.5..1.5]);

unassign('x','y','z'); divf:=diff(force[1], x)+diff(force[2],y)+diff(force[3],z): print(`The divergence of this force is`, divf); totalflux:=int(int(int(divf, z=-sqrt(1-x^2-y^2)..1), y=-sqrt(1-x^2)..sqrt(1-x^2)), x=-1..1): print(`The value of the divergence integral is`, totalflux);

<Text-field style="Heading 1" layout="Heading 1">Selection of Closed Region to Maximize Flux</Text-field>

In this problem, you are to find the dimensions of a particular rectangular box for which the total flux of a force with components LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNiotSSNtaUdGJDYjUSFGJy1GIzYoLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y/LyUpc3RyZXRjaHlHRj8vJSpzeW1tZXRyaWNHRj8vJShsYXJnZW9wR0Y/LyUubW92YWJsZWxpbWl0c0dGPy8lJ2FjY2VudEdGPy8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRk4tRiM2JC1JJW1zdXBHRiQ2JS1GMTYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjtRJ2l0YWxpY0YnLUkjbW5HRiQ2JFEiMkYnRjovJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjotRjc2LVEoJm1pbnVzO0YnRjpGPUZARkJGREZGRkhGSkZMRk8tRiM2KC1GaW42JFEiNEYnRjotRjc2LVExJkludmlzaWJsZVRpbWVzO0YnRjpGPUZARkJGREZGRkhGSi9GTVEmMC4wZW1GJy9GUEZbcEZWRmdvLUYxNiVRInlGJ0ZZRmZuRjpGMEY6LUY3Ni1RIixGJ0Y6Rj0vRkFGZW5GQkZERkZGSEZKRmpvL0ZQUSwwLjMzMzMzMzNlbUYnLUYjNiZGNi1GIzYoLUZpbjYkUSI2RidGOkZnb0ZdcEZnby1GMTYlUSJ6RidGWUZmbkY6RjBGOkZgcC1GIzYmLUZpbjYkUSMxMkYnRjpGZ29GXXFGOkYwRjpGOi8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Y6 is maximized.

<Text-field style="Heading 1" layout="Heading 1">Visualizing the Problem</Text-field>

The force here is the same as in the previous part, but we are now asked to find the coordinates of a rectangular box that will maximize the outward flux. Before we apply the Divergence Theorem, let's plot the figure. Before we apply the Divergence Theorem, let's plot the figure. Remember, you can view the plot from different angles by clicking on the plot and draggin the mouse.

unassign('x','y','z','t','force'); force:=[-x^2-4*x*y, -6*y*z, 12*z]; forceplot:=fieldplot3d(force, x=0..3, y=0..2, z=0..1, labels=["x","y","z"], axes=boxed, grid=[ 6,6,6], arrows=THICK): display(forceplot);

<Text-field style="Heading 1" layout="Heading 1">Evaluating and Maximizing the Divergence Integral </Text-field>

divf:=diff(force[1], x)+diff(force[2],y)+diff(force[3],z): divergenceintegral:=int(int(int(divf, z=0..1), y=0..b), x=0..a);

We will find the values of a and b that maximize this function by computing the two first partials with respect to a and b and equating them both to 0.

ma:=diff(divergenceintegral,a): mb:=diff(divergenceintegral,b): solab:=solve({ma=0, mb=0},{a,b});

Based on these results, what answer will you give? Evaluate the divergence integral at these values for a and b.

solab[1]; eval(divergenceintegral, solab[1]);

solab[2]; eval(divergenceintegral, solab[2]);

solab[3]; eval(divergenceintegral, solab[3]);

solab[4]; eval(divergenceintegral, solab[4]);

You should check this integral out for other values of a and b to convince yourself that you have maximized the divergence integral. Does the fact that this is less than the divergence integral found in Part II mean that you have not maximized the divergence integral as requested?

JSFH

JSFH