Visualizing and Interpreting the Divergence Theorem

Introduction

OBJECTIVE:  Observe that surface integrals are difficult to evaluate, even using Mathematica. 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 Mathematica, 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 Mathematica can adjust the code to solve other problems.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Visualizing the Problem

We select the object formed by a hemisphere of radius 1, topped with a cylinder of radius 1 and height 1, with a circular disk on top on the cylinder. This object is similar to the one to which the method of moments was applied in the previous chapter. The force we choose has components {-(x^2)-4x y,-6 y z, 12z}.

Clear[x,y,z,t,force]
force:={-(x^2)-4x y,-6 y z, 12z}

We will now graph our surface with the force field.

Clear[x,y,z,u,v,t,r,θ]
forceplot=VectorPlot3D[force,{x,-40,40},{y,-40,40},{z,-50,50},VectorPoints->7,VectorScale->{Small,Automatic,#3&}];hemisplot=RevolutionPlot3D[-(Sqrt[1600-(r^2)]),{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
cylpts=Table[{40 Cos[t],40 Sin[t],u},{t,0,2 π,(π/12)},{u,0,40}];
cylplot=ListPointPlot3D[cylpts];
topplot=RevolutionPlot3D[40,{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
flagpolepts=Table[{0,0,t},{t,35,90}];
flagpts=Table[{0,u,v},{u,0,15},{v,70,90}];
fp1=ListPointPlot3D[flagpolepts];
fp2=ListPointPlot3D[flagpts];
Show[forceplot,hemisplot,cylplot,topplot,fp1,fp2,PlotRange->All]

Part II: Evaluating the Divergence Integral

The divergence integral over the volume is easy to set up and calculate. We begin by first loading a package that can perform certain vector computations, including the divergence, curl, and cross product. Be certain to execute the Needs cell first. No output is yielded for Needs.

Needs["VectorAnalysis`"]

Clear[x,y,z]
Print["The divergence of this force is ",divf=Div[force,Cartesian[x,y,z]]]
Print["The value of the divergence integral is ",totalflux=Integrate[divf,{x,-1,1},{y,-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{z,-(Sqrt[1-(x^2)-(y^2)]),1}]]

\[MathematicaIcon] About Mathematica

Part III: Finding Surface Integrals

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.

Hemisphere

Clear[x,y,z,r,θ,ϕ]
x = Sin[θ]Sin[ϕ];
y=Cos[θ]Sin[ϕ];
z=Cos[ϕ];
r={x,y,z};	
rtheta=D[r,θ];
rphi=D[r,ϕ];
cr=Cross[rtheta,rphi]//Simplify;
Print["The integrand for the hemispherical surface integral is ",integrand=force.cr//Simplify]
Print["The value of the hemispherical surface integral is ",hemisflux=Integrate[integrand,{θ,0,2π},{ϕ,π/2,π}], "or ",
nhemisflux=N[hemisflux]]

The equations of the cylinder and the top can be written using cylindrical coordinates, with z and θ the parameters for the cylinder and  r = 1.

Cylinder

Clear[x,y,z]
x=Cos[θ];
y=Sin[θ];
z=z;
r={x,y,z};	
rtheta=D[r,θ];
rz=D[r,z];
cr=Cross[rtheta,rz]//Simplify;
Print["The integrand for the cylindrical surface integral is ",integrand=force.cr//Simplify]
Print["The value of the cylindrical surface integral is ",cylflux=Integrate[integrand,{θ,0,2π},{z,0,1}]," or ",
ncylflux=N[cylflux]]

Top

Polar coordinates r and θ 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.

Clear[x,y,z,r]
x=Cos[θ];
y=Sin[θ];
z=1;
r={x,y,z};	
normal={0,0,1};
Print["The integrand for the surface integral for the top is ",integrand=force.normal//Simplify]
Print["The value of the surface integral for the top is ",topflux=Integrate[R integrand,{θ,0,2π},{R,0,1}]," or ",
ntopflux=N[topflux]]

Part IV: Verifying the Divergence Theorem

Compare your answers to verify the Divergence Theorem.

Print["The total of the surface integrals is ",totalsurfaceflux=hemisflux+cylflux+topflux,"or ",N[totalsurfaceflux]]
If[totalsurfaceflux<=totalflux+.001&&totalsurfaceflux>=totalflux-.001,Print["The Divergence Theorem is validated."],Print["there's a problem"]]

You Try It

Visualize and Evaluate Other Divergence Integrals

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

Clear[x,y,z,u,v,t,r,θ,force]
force:={Cos[10z], Sin[10z],x y} 
forceplot=VectorPlot3D[force,{x,-40,40},{y,-40,40},{z,-50,50},VectorPoints->8,VectorScale->{Small,Automatic,#3&}];hemisplot=RevolutionPlot3D[-(Sqrt[1600-(r^2)]),{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
cylpts=Table[{40 Cos[t],40 Sin[t],u},{t,0,2 π,(π/12)},{u,0,40}];
cylplot=ListPointPlot3D[cylpts];
topplot=RevolutionPlot3D[40,{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
flagpolepts=Table[{0,0,t},{t,35,90}];
flagpts=Table[{0,u,v},{u,0,15},{v,70,90}];
fp1=ListPointPlot3D[flagpolepts];
fp2=ListPointPlot3D[flagpts];
Show[forceplot,hemisplot,cylplot,topplot,fp1,fp2,PlotRange->All]

Now, see if Mathematica can find the divergence integral. if you have already loaded the VectorAnalysis package, there is no need to reload it.

Needs["VectorAnalysis`"]

Clear[x,y,z]
Print["The divergence of this force is ",divf=Div[force,Cartesian[x,y,z]]]
Print["The value of the divergence integral is ",totalflux=Integrate[divf,{x,-1,1},{y,-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{z,-(Sqrt[1-(x^2)-(y^2)]),1}]]

Selection of Closed Region to Maximize Flux

In this problem, you are to find the dimensions of a particular rectangular box for which the total flux of a force with components {-(x^2)-4x y,-6 y z, 12z}  is maximized.

Visualizing the Problem

If you have not already loaded the packages used below, be sure to do that. The force here is the same as in Part I, 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.

Clear[x,y,z,t,force]
force:={-(x^2)-4 x y,-6 y z,12 z}
forceplot=VectorPlot3D[force,{x,-40,40},{y,-40,40},{z,-50,50},VectorPoints->8,VectorScale->{Small,Automatic,#2&},AxesLabel->{x,y,z}]

Evaluating and Maximizing the Divergence Integral

divf=Div[force,Cartesian[x,y,z]];
divergenceintegral=Integrate[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=D[divergenceintegral,a];
mb=D[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.

divergenceintegral/.solab

You should check this divergence integral out for other values of a and b and see if the your result is validated. Also, how is this value for the divergence integral related to the value you arrived at in Part II for this same force? Does the fact that that answer is larger contradict your having found the maximum here?

\[MathematicaIcon] About Mathematica

In computing multiple integrals, Mathematica starts with the variable and bounds on the extreme right of the Integrate command and then continues the evaluation from right to left.
Go back.